A spectral Galerkin method for the the coupled Orr-Sommerfeld 
and induction equations for free-surface MHD 

Dimitrios Giannakis a> * , Paul F. Fischer b , Robert Rosner a,b ' c 

a Department of Physics, University of Chicago, Chicago, IL 60637, USA 
b Argonne National Laboratory, Argonne, IL 60439, USA 
c Department of Astronomy and Astrophysics, 
University of Chicago, Chicago, IL 60637, USA 



Abstract 

We develop and test spectral Galerkin schemes to solve the coupled Orr-Sommerfeld (OS) and induction equations 
for parallel, incompressible MHD in free-surface and fixed-boundary geometries. The schemes' discrete bases consist 
of Legendre internal shape functions, supplemented with nodal shape functions for the weak imposition of the 
stress and insulating boundary conditions. The orthogonality properties of the basis polynomials solve the matrix- 
coefficient growth problem, and eigenvalue-eigenfunction pairs can be computed stably at spectral orders at least 
as large as p — 3,000 with p-independent roundoff error. Accuracy is limited instead by roundoff sensitivity due 
to non-normality of the stability operators at large hydrodynamic and/or magnetic Reynolds numbers (Re, Rm > 
4 x f0 4 ). In problems with Hartmann velocity and magnetic-field profiles we employ suitable Gauss quadrature rules 
to evaluate the associated exponentially weighted sesquilinear forms without error. An alternative approach, which 
involves approximating the forms by means of Legendre-Gauss-Lobatto (LGL) quadrature at the 2p — 1 precision 
level, is found to yield equal eigenvalues within roundoff error. As a consistency check, we compare modal growth 
rates to energy growth rates in nonlinear simulations and record relative discrepancy smaller than 10 -5 for the least 
stable mode in free-surface flow at Re = 3 x I0 4 . Moreover, we confirm that the computed normal modes satisfy 
an energy conservation law for free-surface MHD with error smaller than 10 -6 . The critical Reynolds number in 
free-surface MHD is found to be sensitive to the magnetic Prandtl number Pm, even at the Pm — O(10~ 5 ) regime 
of liquid metals. 

Key words: Eigenvalue problems, spectral Galerkin method, hydrodynamic stability, Orr-Sommerfeld equations, free-surface 
MHD 
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1. Introduction 



The Orr-Sommerfeld (OS) and induction equations, Eqs. (2.1) below, govern the linear stability of tempo- 
ral normal modes in incompressible, parallel magnetohydrodynamics (MHD). These equations have mainly 
been applied to study the stability of flows with fixed domain boundaries in the presence of an external mag- 
netic field ([1] and references therein). However, linear-stability analyses of free-surface flows have received 
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comparatively little attention. Here the OS and induction equations, in conjunction with the kinematic 
boundary condition at the free surface (2.5), pose a coupled eigenvalue problem which must be solved for 
the complex growth rate 7, the velocity and magnetic-field eigcnfunctions, respectively u and b, as well as 
the free-surface oscillatory amplitude a. 

Free-surface MHD arises in a variety of contexts, including cooling of fusion reactor walls by liquid-metal 
blankets [2], liquid-metal forced flow targets [3], and surface models of compact astrophysical objects [4]. In 
these and other cases of interest, hydrodynamic Reynolds numbers are large (Re > 10 4 ), and the flow takes 
place in the presence of a strong background magnetic field (H > 10 2 , where H is the Hartmann number). 
All terrestrial fluids have small magnetic Prandtl numbers (e.g. for laboratory liquid metals Pm < 10~ 5 ), 
suggesting that the magnetic field is well in the diffusive regime. On the other hand, Pm = 0(1) flows have 
been conjectured to play a role in certain astrophysical accretion phenomena [5]. 

The main objective of the present work is to develop accurate and efficient spectral Galerkin schemes 
for linear-stability analyses of free-surface and fixed-boundary MHD. Our schemes build on the Galerkin 
method for plane Poiseuillc flow by Kirchner [6], and Mclcnk, Kirchncr, and Schwab [7], hereafter collectively 
referred to as KMS. A companion article [8] discusses the operating physics in low-Pm, problems. A future 
objective is to test our linear models against wave-dispersion and critical-parameter data from a free-surface 
MHD experiment at Princeton Plasma Physics Laboratory (PPPL) by Ji and coworkers [9,10]. 

1.1. Background 

Since the pioneering work of Orszag [11] in 1971, spectral methods have emerged as a powerful tool to 
solve hydrodynamic-stability problems. Orszag applied a Chebyshev tau technique to transform the OS 
equation for plane Poiseuillc flow to a matrix generalized eigenproblem Ku = jMu, which he solved at 
Reynolds numbers of order 10 4 using the QR algorithm. The superior performance of the Chebyshev tau 
method compared to existing finite-difference and spectral schemes led to its application to a diverse range 
of stability problems (e.g. [12]). However, despite the widespread use of spectral techniques in flows with 
fixed domain boundaries, most numerical stability analyses of free-surface problems to date are based on 
finite-difference methods. Among these are the studies of gravity and shear-driven flows by Dc Bruin [13], 
and Smith and Davis [14]. To our knowledge, the only related work in the literature employing spectral 
techniques is contained in the PhD thesis by Ho [15], where the OS equation for a vertically falling film is 
solved at small Reynolds numbers (Re < 10). 

In MHD, numerical investigations on the stability of modified plane Poiseuille flow subject to a transverse 
magnetic field, also known as Hartmann flow, begin in 1973 with the work of Potter and Kutchey [16], who 
used a Runge-Kutta technique to solve the coupled OS and induction equations at small Hartmann numbers 
(H < 6). Lingwood and Alboussiere [17] also employed a Runge-Kutta method to study the stability of an 
unbounded Hartmann layer. An early application of spectral methods was performed by Dahlburg, Zang, 
and Montgomery [18] in 1983, who adopted Orszag's scheme to investigate the stability of a magnetostatic 
quasiequilibrium (i.e. a state where the fluid is at rest but a slowly varying background magnetic field 
is present). A Chebyshev tau method for plane Poiseuille and plane Couette flows in the presence of a 
transverse magnetic field was later developed by Takashima [19,20]. Takashima's calculations extend to high 
Reynolds and Hartmann numbers (Re ~ 10 7 , H ~ 10 3 ) and over a range of magnetic Prandtl numbers up 
to Pm = 0.1. In addition, he considers the limiting case of vanishing magnetic Prandtl number, where the 
OS and induction equations are replaced by a single equation (2.2). However, his analysis does not take into 
account modes other than the least stable one (cf. [6,11,12]). 

A major challenge in hydrodynamic-stability problems at high Reynolds numbers is the existence of thin 
boundary layers, whose thickness scales as (aRe)~ x / 2 for a normal mode of wavenumber a [21], which 
requires the use of large spectral orders p to achieve convergence. Specifically, Melenk et al. [7] have shown 
that a necessary condition for accurate results is that the ratio Re/p 2 is small, implying that for problems 
of interest the required p can be in the thousands. At such high spectral orders the Chebyshev tau method 
can be problematic, since it gives rise to stiffness and mass matrices, respectively K and M, that are (i) 
densely populated (the storage and computation cost therefore scale as p 2 and p 3 , respectively), and (ii) ill 
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conditioned (the matrix elements associated with a fourth-order differential operator, such as the OS one, 
grow as p 7 ). One way to alleviate the matrix-coefficient growth problem is to pass to a streamfunction- 
vorticity formulation [22], or, more generally, apply the D 2 method proposed by Dongarra, Straughan, and 
Walker [12]. Here one achieves a p 3 coefficient scaling by casting the OS equation into two coupled second- 
order equations, but at the expense of doubling the problem size. 

Another drawback of the tau method is the occurrence of 'spurious' eigenvalues, i.e. eigenvalues with large 
magnitude (e.g. O(10 17 ) [23]), and real-part oscillating between positive and negative values as p is varied. 
These numerical eigenvalues are not at all related to the spectrum of the OS operator, and in order to avoid 
drawing erroneous conclusions (e.g. deciding that a flow is unstable when the unstable mode is spurious), 
the practitioner must either detect them and ignore them in the analysis (the non-spurious modes are 
computed correctly), or eliminate them by a suitable modification of the method (e.g. [12,22]). Their origin 
has been elucidated by Dawkins, Dunbar, and Douglas [24], who found that the large spurious eigenvalues 
in Chebyshev tau schemes are perturbations of infinite eigenvalues in nearby Legendre tau discretizations. 

Recently, KMS have developed a spectral Galerkin method that addresses some of the aforementioned 
shortcomings. Central to their scheme is the use of the so-called compact combinations of Legendre poly- 
nomials [25,26], or hierarchical shape functions [27], as a basis of the Sobolev space H% (the trial and test 
space for velocity eigenfunctions). The resulting orthogonality properties solve the matrix-coefficient growth 
problem, and no reduction in the differential-equation order is required (in fact, the condition number of 
K has been found to be independent of p > 100 [7]). Moreover, the stiffness and mass matrices are sparse, 
provided that the basic velocity profile is polynomial. In that case, memory requirements scale as p, and 
iterative solvers can be used to compute the eigenvalues and eigenvectors efficiently. A further attractive 
feature of the method, which appears to be connected to the non-singularity of M, is that it gives no rise 
to spurious eigenvalues. 

An additional, and perhaps more fundamental, challenge is related to the non-normality of hydrodynamic- 
stability operators, which becomes especially prominent at large Reynolds numbers. In that regime, even 
though the eigenfunctions may form a complete set (as has been proved for bounded domains [28]), they are 
nearly linearly dependent. A key physical effect of the eigenfunctions' non-orthogonality is large transient 
growth of asymptotically stable perturbations, which suggests that eigenvalue analysis is of little physical 
significance [29] . An alternative method that aims to capture the effects of transient growth is pseudospec- 
tra [30], but this will not be pursued here. We remark, however, that comparisons between spectra and 
pseudospectra are common in pseudospectral analyses, and stable and efficient schemes for eigenvalue com- 
putations are desirable even in that context. 

At the numerical level, non-normality is associated with high sensitivity of the spectrum to roundoff 
errors [31]. This effect was noted by Orszag himself [11], who observed significant changes in the computed 
eigenvalues by artificially reducing numerical precision from 10~ 14 to 10~ 8 . The eigenmodes that are most 
sensitive to perturbations of the OS operator and its matrix discrete analog are those lying close to the 
intersection point between the A, P, and S eigenvalue branches on the complex plane (see [32] for a description 
of the nomenclature) . Reddy, Schmid, and Henningson [33] have observed that in plane Poiseuille flow at 
Re <~ 10 4 perturbations of order 10~ 6 suffice to produce 0(1) changes in the eigenvalues near the branch 
point. Moreover, they found that roundoff sensitivity increases exponentially with the Reynolds number. 
Qualitatively, this type of growth is attributed to the existence of solutions of the OS equation that satisfy 
the boundary conditions to within an exponentially small error. In consequence, double-precision arithmetic 
(typically 15 digits) rapidly becomes inadequate, and for Re > 4 x 10 4 one obtains a diamond shaped pattern 
of numerical eigenvalues instead of a well resolved branch point (e.g. Fig. 4 in [12] and Figs. 14-16 below). 
Dongarra et al. [12] have postulated that alleviation of this spectral instability requires the use of extended- 
precision arithmetic, and cannot be removed by improving the conditioning of the numerical scheme (e.g. 
employing a D method instead of D 2 one). Melenk et al. [7] note that their Galerkin method accurately 
resolves the eigenvalue branch point at Re = 2.7 x 10 4 using 64-bit arithmetic, when the Chebyshev tau 
method already produces the diamond-shaped pattern. However, even a moderate Reynolds-number increase 
(e.g. Re = Ax 10 4 in Fig. 14) results to the appearance of the pattern, despite the Galerkin scheme's superior 
conditioning. It therefore appears that, at least in these examples, the decisive factor in roundoff sensitivity 
is the non-normality of the OS operator rather than the details of the discretization scheme. 
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1.2. Plan of the Present Work 



The principal contribution of this article is twofold. First, we generalize the spectral Galerkin method of 
KMS for plane Poiseuille flow to free-surface and fixed-boundary MHD. Second, we present a number of test 
calculations aiming to assess our schemes' numerical performance, as well as to provide benchmark data. 
The calculations have been performed using a Matlab code, available upon request from the corresponding 
author. 

As already stated, central to the stability and efficiency of the KMS scheme is the use of suitable linear 
combinations of Legendre polynomials as a basis of Hq. In the sequel, we employ similar constructions to 
treat the free-surface MHD problem. Here the velocity field obeys stress conditions at the free surface, which 
we enforce weakly (naturally) by supplementing the basis with nodal shape functions [27]. Pertaining to the 
magnetic field, we assume throughout that the domain boundaries are electrically insulating, from which 
it follows that it obeys boundary conditions of Robin type, with extra contributions from the free-surface 
oscillation amplitude [8]. We enforce naturally these boundary conditions as well, discretizing the magnetic- 
field by means of the internal and nodal shape functions for H 1 . As we demonstrate in §5.2.1 and §5.2.2, 
our choice of bases gives rise (and is essential) to a major advantage of our schemes, namely that roundoff 
error is independent of the spectral order p. 

In problems with polynomial steady-state profiles the stiffness and mass matrices are sparse, and closed- 
form expressions exist for their evaluation (see Appendix A). On the other hand, in Hartmann flow K 
becomes full and must be computed numerically, since the discretization procedure introduces inner products 
of Legendre polynomials with exponential weight functions. We evaluate the required integrals stably and 
without error by means of the Gauss quadrature rules developed by Mach [34], who has studied a class of 
orthogonal polynomials with exponential weight functions on a finite interval. Following the standard practice 
in finite-element and spectral-element methods [35,36], we also consider a method where the problem's 
weighted sesquilinear forms are replaced by approximate ones derived from Legendre-Gauss-Lobatto (LGL) 
quadrature rules. At an operational level, the latter approach has the advantage of being sufficiently general 
to treat arbitrary analytic profiles. However, it introduces quadrature errors, and one has to ensure that the 
stability and convergence of the scheme are not affected. As shown by Bancrjee and Osborn [37], in finite- 
element schemes for elliptical eigenvalue problems that is indeed the case, provided that the approximated 
eigenfunctions are smooth and the quadrature rule is exact for polynomial integrands of degree 2p — 1. 
To our knowledge, however, no such result exists in the literature for the OS eigenproblcms we study 
here, and is therefore not clear what (if any) quadrature precision would suffice. Even though we make 
no attempt to parallel Banerjee and Osborn's work, we nevertheless find that eigenvalues computed using 
approximate quadrature at the 2p — 1 precision level converge, modulo roundoff error, to the same value as 
the corresponding ones from the exact-quadrature method. 

One of the advantages of the spectral Galerkin method is its flexibility. Our scheme for free-surface MHD 
can be straightforwardly adapted to treat MHD problems with fixed domain boundaries, problems in the 
limit of vanishing magnetic Prandtl number, as well as non-MHD problems. In §5.1 we describe the basic 
properties of the eigenvalue spectra of these problems, leaving a discussion of the physical implications 
to Ref. [8]. We also present a series of critical-parameter calculations (see §5.4), confirming that results 
obtained via the fixed-boundary variants of our schemes are in close agreement with the corresponding ones 
by Takashima [19]. In free-surface problems, when Pm is increased from 1CP 8 to 10~ 4 the critical Reynolds 
number is seen to drop by a factor of five, while the corresponding relative variation in fixed-boundary 
problems is less than 0.003. Due to the limited availability of eigenvalue data for free-surface flow (cf. fixed- 
boundary problems [6,7,12,19]), we were not able to directly compare our free-surface schemes to existing 
ones in the literature. Instead, we have carried out two other types of consistency checks (see §5.3), one of 
which is based on energy-conservation laws in free-surface MHD, whereas the second involves growth-rate 
comparisons with fully nonlinear simulations. 

A numerical caveat concerns the aforementioned roundoff sensitivity at high Reynolds numbers. In §5.2.3 
we observe that as Re grows our schemes experience the spectral instability that has been widely encountered 
in Poiseuille flow [6,7,11,12,33]. Most likely, this issue is caused by the physical parameters of the problem, 
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rather than the properties of the discretization scheme, and can only be addressed by increasing the numerical 
precision. Unfortunately, since the latter option is (as of January 2008) not natively supported in Matlab, we 
merely acknowledge the existence of the instability, and work throughout in double-precision arithmetic. We 
remark, however, that only the eigenvalues near the branch-intersection points are affected. In particular, 
eigenvalues and cigcnfunctions at the top end of the spectrum can be accurately computed at Reynolds 
numbers at least as high as 10 . We also wish to note that the emphasis of our work is towards the numerical, 
rather than analytical, side, and even though techniques to study the stability and convergence of Galcrkin 
methods for eigenvalue problems are well established in the literature ([38] and references therein), we do 
not pursue them here. 

The plan of this paper is as follows. In §2 we specify the governing equations and boundary conditions of 
our models. In §3 we develop their weak formulation. The associated Galcrkin discretizations are described in 
§4. We present our numerical results in §5, and conclude in §6. Appendix A contains closed form expressions 
for the matrix representations of the sesquilinear forms used in the main text. Although some of these can 
also be found in [6] , we opted to reproduce them here because that paper contains a number of typographical 
errors. Finally, in Appendix B we have collected tables of eigenvalues for the problems examined in §5.1. 

2. Problem Description 

2.1. Governing Equations 

Using x and z to denote the streamwise and flow-normal coordinates, and D to denote differentiation with 
respect to z, we consider the coupled OS and induction equations, 

Re-^D 2 -a 2 fu - (7 + iaU){V 2 -a 2 )u + ia(D 2 U)u + (iaB x + B z D)(D 2 -a 2 )b - k*(D 2 B x )b = 0, (2.1a) 

and 

Rm- l {T> 2 -a 2 )b - (7 + iaU)b + (iaB x + B z D)u = 0, (2.1b) 

defined over an interval Q = (z\, Z2) € K. Here u € C 4 (J?) and b e C 3 (H) are respectively the velocity and 
magnetic-field eigenfunctions corresponding to the eigenvalue 7 € C. Also, a > is the wavenumber, and 
Re > and Rm > are the hydrodynamic and magnetic Reynolds numbers. The functions U <G C 2 (D) 
and B x <G C 2 (i7) are the steady-state velocity and streamwise magnetic field. The flow-normal, steady-state 
magnetic field B z is constant, since (B X ,B Z ), where (-,-) stands for (x,z) vector components, is divergence 
free and streamwise invariant. The two-dimensional velocity and magnetic fields associated to u and b are 
given by Re((i D u/a, u)e lax+lt ) and Re((iDb/a,b)e iax+ ^). 

A physical derivation of (2.1), as well as (2.2) and the boundary conditions (2.3)-(2.10) ahead, can be 
found in Refs. [8,19]. Here we note that the magnetic-field variables b, B x , and B z have been rendered to 
non-dimensional form using the characteristic magnetic-field Bo := (popY^Uo, where [iq, p, and Uo are the 
permeability of free space, the fluid density, and the characteristic velocity, respectively. With this choice of 
magnetic- field scale, u and b arc naturally additive. Another option (employed e.g. by Takashima [19]) is to 
set Bq = B' , where B' is the typical steady-state magnetic field. The resulting non-dimensional eigenfunction 
b' is related to the one used here according to b' = Ab, where A := (iiqpY^Uq/ B' is the Alfvcn number. 
We also remark that we have adopted the eigenvalue convention used by Ho [15], under which Re(7) =: r 
corresponds to the modal growth rate (i.e. a mode is unstable if r > 0), while C := — Im(7)/a is the phase 
velocity. The complex phase velocity c = i^y/a, where Re(c) = C and Im(c)a = r, is frequently employed in 
the literature (e.g. [6,7,12-14,19]) in place of 7. 

Let Pm := Rm/ Re denote the magnetic Prandtl number; the ratio of magnetic to viscous diffusivity (see 
e.g. [39] for an overview of the dimensionless parameters in MHD). A limit of physical interest, hereafter 
referred to as the inductionless limit [1], is Pm — > with B x independent of z, and the Hartmann numbers 
H x := B x RePm 1/2 and H z := B z RePm 1/2 , measuring the square root of the ratio of Lorentz to viscous 
forces, non-negligible. This situation corresponds to a fluid of sufficiently high magnetic diffusivity so that 
magnetic-field perturbations are small (| \b\ | <C \ \u\ \ in some suitable norm), but Lorentz forces due to currents 
induced by the perturbed fluid motions u within the steady-state field (B X ,B Z ) are nonetheless present. It 
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Fig. 1. Geometry of channel (left) and film (right) problems. U(z) and B(z) are the steady-state velocity and induced magnet- 
ic-field profiles, respectively (see §2.3). 

is then customary to make the approximation (D 2 — a 2 )b = —Rm(\aB x + B z D)u [40] and replace (2.1) by 
the single equation 

(D 2 -c?fu - (iaH x + H z Dfu - Re(<y + iaU){B 2 -a 2 )u + iaRe(D 2 U)u = 0. (2.2) 
2.2. Boundary Conditions 

We study two types of problems, which we refer to as channel and film problems according to their geomet- 
rical configuration (see Fig. 1 for an illustration). Within each category we further distinguish among MHD 
problems, their counterparts in the inductionlcss limit, and non-MHD problems, where all electromagnetic 
effects are neglected. 

In channel problems the flow takes place between two fixed, parallel plates. As is customary, we make the 
domain choice fi = fi c := (—1, 1), and enforce the no-slip boundary conditions 

u(±l) = Du(±l) = 0. (2.3) 

Moreover, we assume that the plates and the region exterior to the flow are perfect insulators, which leads 
to the Robin boundary conditions 

D6(±l)±a6(±l) = (2.4) 

for the magnetic field. 

In film problems we set fi = fit := (—1,0), and consider that the domain boundary z 2 = corresponds 
to a free surface, whose oscillation amplitude a <G C obeys the kinematic boundary condition 

u(0) - (7 + iaU(0))a = 0. (2.5) 

We assume that the free surface is subject to surface tension and gravity (see e.g. [41] for a discussion of 
free-surface dynamics). The non-dimensional stress due to surface tension is given by act 2 / (OhRe) 2 , where 
the Ohncsorge number Oh measures the ratio of viscous to capillary velocity scales. Oh is related to the 
Weber number We (the ratio of surface tension to inertial stresses) via Oh = 1/ (ReWe 1 ^ 2 ). Moreover, 
we express the z component of the gravitational force as — (PgRe)~ 2 , where Pg is a parameter that we 
refer to as the gravitational Prandtl number. Pg is equal to the ratio between the kinematic viscosity v 
and the diffusion constant (t/cos(0)Z 3 ) 1 / 2 , formed by the flow-normal gravitational-field strength, gcos(9), 
and the characteristic length scale I (cf. the magnetic Prandtl number Pm = v/rj, where r\ is the magnetic 
diffusivity) , and is related to the Froude number Fr (the ratio of convective to gravity velocity scales) 
according to Pg = Fr/ Re [8]. Our use of the parameters Oh and Pg, rather than the more familiar We and 
Fr, is motivated by the fact that they do not depend on the characteristic flow velocity, and are thus likely 
to remain fixed in situations where a single working fluid is driven at different flow speeds. Typical values 
for a laboratory liquid-metal film of thickness ~ 1 cm are Oh ~ 10~ 4 and Pg ~ 10 -4 (e.g. [42-44]). 
Balancing the forces acting on the free surface leads to the normal-stress condition 

(((D 2 -3a 2 ) D -i?e( 7 + iaU) D +i<xRe(D U))u)\ z=0 + Re(B z (D 2 -a 2 ) - ia(D B x ))b\ z=0 

a " {pghfe + ~WWe + ** BM ° B ' M 2[aBU ^) a = °' ^ 
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and the shear-stress condition 

D 2 w(0) + a 2 w(0)-k*D 2 [/(0)a = 0. (2.7) 
The no-slip boundary conditions are again 

u(-l) = Du(-l) = 0, (2.8) 

but the insulating boundary conditions 

D6(-l) -ab(-l) = D6(0) + ab(0) - iaD B x (0)a = (2.9) 

now involve the free-surface oscillation amplitude (cf. (2.4)). In the inductionlcss limit, the boundary con- 
ditions for b are not required. Furthermore, (2.6) reduces to 

((D 2 -3a 2 ) D -#e(7 + iaU) D +iaRe(B U) - H z (iaH x + H z B))u\ z=0 

-a 2 ( — I — + — ^ 2iaBU(0)] a = 0. (2.10) 

\Pg 2 Re Oh 2 Re J 

To summarize, we refer to all problems involving a free surface as film problems and those that take 
place within fixed boundaries as channel problems. Within the film category, we call film MHD problems 
those governed by (2.1), subject to the boundary conditions (2.5)-(2.9). These are to be distinguished from 
inductionless film problems, where the coupled OS and induction equations are replaced by (2.2), and the 
boundary conditions arc (2.5), (2.7), (2.8), and (2.10). Similarly, we differentiate between channel MHD 
problems, specified by (2.1), (2.3), and (2.4), and their inductionless variants, where the differential equation 
and boundary conditions are respectively (2.2) and (2.3). Finally, for what we refer to as non-MHD film 
problems and non-MHD channel problems we set H x = H z = in (2.2) and (2.10). We mention in passing 
that one can treat in a similar manner 'jet' problems, where free-surface boundary conditions are enforced 
at z — ±1, although problems of this type will not be considered here. 

2.3. Steady-State Configuration 

In what follows we consider the magnetic-field configuration 

(B x (z),B z ) = {A-\A- Z l ) + (A-'RmBiz), 0), (2.11) 

where (A" 1 , A J 1 ) is a uniform, externally imposed magnetic field, quantified in terms of the streamwise 
and flow- normal Alfvcn numbers A x and A z , and B <G C 2 (i?) is a function representing the magnetic field 
induced by the fluid motion U (z) within the background field (B is equal to the corresponding function B 
in [19]). For the test calculations presented in §5 we employ the Hartmann profiles [1] 

U(z) = (cosh(H z ) - cosh(H z z))/X, H z B(z) = {smh(H z z) - smh(H z )z)/X, (2.12) 

where H z = (ReRm) 1 / 2 A~ x , X = cosh(iJ z ) — 1, and z € [—1, 1]. Note that the expressions above are valid 
for both channel and film problems. In the latter case, one restricts z to the interval Df to obtain 'half of 
the corresponding channel profile. A further useful quantity is the mean velocity, 

(10 := f * dz = (cosh(H z ) - sinh(H z )/H z )/X, (2.13) 

which grows monotonically from 2/3 to 1 as H z increases from zero to infinity. The steady-state configuration 
described by (2.11) and (2.12) is a solution of the unperturbed Navier-Stokes and induction equations [8]. 
In the limit H z — > we have, 

U(z) = l-z 2 , B(z) = -z(l-z 2 )/3, (2.14) 

indicating that the velocity profile reduces to the usual Poiscuillc one. Even though B is nonzero in the 
limit, the streamwise induced magnetic field A~ x RmB = Pm}^ 2 H z B vanishes. For H z > the velocity 
and magnetic-field profiles develop exponential tails of thickness 1 /H z , where the vorticity and current are 
concentrated. These so-called Hartmann layers form near the no-slip walls, as shown in Fig. 2. 
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Fig. 2. Steady-state velocity U (left) and magnetic field B (right) for Hartmann flow (2.12) at H z = 0, 10, 20 

2.4. Energy Balance 

In §5.1 and §5.3 ahead we shall make use of energy conservation laws for the normal modes, which follow 
from the linearized Navier-Stokes and induction equations governing the evolution of linear perturbations 
in MHD. Leaving the details of the derivation to [8], to each (u, b, a) satisfying (2.1) and the boundary 
conditions (2.5)-(2.9) we assign an energy E := E u + E h + E a , consisting of kinetic, magnetic, and surface 
contributions 



E u := J dz (\Du(z)\ 2 + a 2 \u(z)\ 2 ), 



(2.15a) 



E b :-- 



dz (\Db(z)\ 2 + a 2 \b(z)\ 2 ) + 2a(\b(-l)\ 2 + \b(0)\ 2 ), 



(2.15b) 

E a := a 2 (Fr- 2 + B x (0) B B x (0) + Wea 2 ) \a\ 2 . (2.15c) 

Here the kinetic energy E u is (up to a proportionality constant) the energy norm of the 2D velocity field 
associated with the velocity eigenfunction u, while the magnetic energy contains, in addition to the energy 
norm of the magnetic field within the fluid, boundary terms representing the energy of the field penetrating 
through the insulating boundaries. The surface energy consists of potential energy due to gravitational and 
magnetic stresses, plus a contribution from surface tension. In inductionless problems the modal energy is 
E = E u + E a , where E u is given by (2.15a), and E a follows from (2.15c) with B x set to zero. 
Aside from E, to each (u, b, a) correspond power-transfer terms 

(2.16a) 



r, 



dz (DU{z))Im(u{z)*Du{z)), 

E J-i 

-| J° dz (BU(z))lm(b(z)*Bb(z)), 

— / dz (DB x (z))lm(u(z)*T)b(z)-b(z)*Du(z)), 
E J -i 

r„:=--^- J dz (\D 2 u(z)\ 2 -2a 2 Rc(u(z)*D 2 u(z)) + a 4 \u(z)\ 2 ), 

l — [ dz (| D 2 b(z)\ 2 - 2a 2 Rc(b(z)* D 2 b(z)) + a A \b(z)\ 2 ), 
Rm J_! 

(DI7(0))Im(Du(0)o*), 



r — 

n.i/ • — 



ERm 

a 

ERe 
a 



(2.16b) 

(2.16c) 

(2.16d) 

(2.16e) 
(2.16f) 



ar > - ERm 



(D J5 X (0)) (lm((D 2 6(0) - a 2 b(0))a*) + B z Im(Du(0)o*) + aB x (0) Re(u(0)a*)) , (2.16g) 
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each of which has a physical interpretation. Jr and Im are the Reynolds and Maxwell stresses, i.e. the 
power transferred from the steady-state velocity field U to the velocity and magnetic-field perturbations. 
Ij is the so-called current interaction; the power transfer from the steady-state current (represented by 
the T)B X term in (2.16c)) to the perturbations. The non-positive quantities F v and r v are respectively the 
viscous and resistive dissipation (the subscripts v and r\ stand for the viscous and magnetic diffusivities) . 
Finally, the surface terms T av and F ar] represent the power transferred to the free surface by viscous and 
electromagnetic forces, respectively. It can be shown that the sum of the terms in (2.16) is equal to the 
modal growth rate. That is, the real part F of the eigenvalue 7 corresponding to (u, b, a) is expressible as 

r = Re( 7 ) = r R + r M + r r + r v + r v + r au + r ar) . (2.17) 

Similar energy-balance relations can be derived for channel and inductionless problems, but we do not require 
them here. 

3. Weak Formulation 

We now cast the eigenvalue problems specified in §2 into weak (variational) form, suitable for the Galcrkin 
schemes developed in §4. With a slight abuse of notation we use the symbol Q to denote the domain of 
both film and channel problems, where it is understood that fl stands for cither i?f (film problems) or fl c 
(channel problems), depending on the context. Also, we collectively denote the vector spaces of admissible 
solutions for the velocity and magnetic- field eigenfunctions by V u and Vb, respectively, even though different 
versions of these spaces will be constructed for film and channel problems. In what follows, we describe the 
procedure of obtaining the weak formulation of film MHD problems. Channel MHD problems, as well as 
the inductionless variants of film and channel problems, can be treated in an analogous manner, and, in the 
interests of brevity, we shall merely state the results. 

Given an interval £2 = (z\, Z2) € K, we denote by L 2 {D) the Hilbert space of square- integrable complex- 
valued functions on J?, equipped with the inner product {v\,V2)o,n '■— dz V\{z)v2{z)* and induced norm 
IMIo.f2 '■— i v > v )o,n- We then introduce as usual (e.g. [45]) the Sobolev spaces H k (f2), k <G N, consisting of 
elements v € L 2 (f2), whose weak derivatives D d v for \d\ < k are also in L 2 {f2). Moreover, H^(f2) are the 
closures in H k (f2) of C§°(i2), the space of smooth, compactly supported functions on Q. The associated 
semi-norms and norms are given by \v\l n := || D fc v\^ n and |Mlfc,f2 := ELo Ml: where | • \ k ,n and || • \\ kt a 
are equivalent norms on Hq(Q). Using the symbol to denote embedding, it is a consequence of the 
Sobolev embedding theorem that H 2 ([2) ^> C 1 (i7) [45]. That is, each v e H 2 {Q) is equal to a unique 
function v € C 1 (/5), except on a measure-zero subset of Q. This allows us to define the boundary- value 
maps S^ : H 2 {Q) C for i e {1, 2} and j € {0, 1}, where S^(w) = D- 7 (w(z i )). We then construct the space 

H 2 {Q) := {v e H 2 {Q)- S?(w) - S\(v) - 0}, (3.1) 

which will be used as trial and test space of velocity eigenfunctions in film problems. Using the embedding 
H X {Q) ^> C°(f2), we also introduce the boundary-value maps S°(u) = v(zi) for H 1 (Q), where now v G 
H 1 (n) and v is its image in C°(£2). The latter two maps will be used to (weakly) enforce the insulating 
boundary conditions obeyed by the magnetic field. 

In the strong (classical) formulation of film MHD problems we express Eqs. (2.1) and (2.5) in the form 

JC(u,b,a) = 'yM(u,b,a), (3.2) 

where IC and M. are matrix differential operators. These so-called 'stiffness' and 'mass' operators, respectively 
with domain V K = C 4 (/5) x C 2 {Q) x C and V M = C 2 {f2) x C 1 (J?) x C D V Kl are given by 





(ic 


JCub 










' -Re(B 2 


-a 2 ) o\ 




(u\ 


K{u, b, a) = 


K-bu 


JCbb 






b 


, M(u, b, a) = 





Rm 




b 




cO 


- 


iaU(0) J 








v° 


1) 




[a) 
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where 

ICuu = -(D 2 -a 2 ) 2 + iaRe(U(D 2 -a 2 ) - (D 2 U)) 7 JC bb = D 2 -a 2 - iaRmU, (3.4a) 
JC ub = -Re(iaB x + B Z ~D)(D 2 -a 2 ) + iaRe(D 2 B x ), K, bu = Rm(iaB x + B z D). (3.4b) 

Note that in (3.2) we have multiplied the OS equation (2.1a) by — 1. This is a conventional manipulation, 
with no influence on the scheme's numerical behavior, made in order to obtain a positive-definite mass 
form in (3.5) below. The strong version of the problem may then be stated as follows: Find 7 e C and 
(u, 6, a) e T>ic \ {(0, 0, 0)}, such that the governing equations (3.2), and the boundary conditions (2.5)-(2.9) 
are satisfied. 

In order to pass from the strong to the weak (variational) formulation, one begins by identifying the 
spaces of admissible solutions V u and V b for the velocity and magnetic-field cigenfunctions, respectively. 
In film problems we set V u = Hf(Q) and V b — i? 1 (/7), so that the no-slip boundary conditions (2.8) are 
enforced strongly, whereas the stress and insulating boundary conditions, (2.6), (2.7), and (2.9), must be 
imposed in a natural (weak) sense (e.g. [15,46]). Taking the free-surface amplitude into account, the full 
solution space is therefore V = V u x V b x C, which we equip with the direct-sum inner product (vi, V2)v,n = 
{ui, u 2 )o,n + (61, b 2 )o,Q + 0102*, where Uj e V u , bj <G V b , aj G C, and Vj = (uj,bj,aj) G V for j e {1, 2}. 

We now proceed to construct scsquilinear forms K and M associated to JC and M. , respectively. Introducing 
a test element v = (u, 6, a) € V, we form the (•, -)v,n inner product of (3.2) with v, namely (IC(u, b, a),v)v,n = 
j(M(u, b, a), v)v,n- Upon integration by parts this leads to 

K(v,v) =-yM(v,v), (3.5) 

where now v = (u,b,a) € V D V{JC). Also, K : V x V C and M : V x V 1— > C are sesquilinear forms 
associated with the mass and stiffness operators K. and M, respectively. We make the decompositions 

K(v, v) = K uu (u, u) + K ub (b, u) + K ua (a, u) + K bu (u, b) + K bb (b, b) + K ba (b, a) 

+ K au (u,a) + K aa (a, a), (3.6a) 
M(v,v) = M uu (u,u) + M bb (b,b) + M aa (a,a), (3.6b) 

which consist of the following objects: In (3.6a), K uu : V u x V u 1— > C, K bb : V b x V b 1— ► C, and K aa :CxC^C 
are sesquilinear forms given by 

K uu (u,u) = K [ °l(u,u) + K^(«,«) + K{S(«,u), (3.7a) 

K bb (b,b) = K [ $(b,b) + K [ b u b ] (b,b) + K™(6,6), (3.7b) 
K aa (a,a) = -\aU(0)aa*, (3.7c) 

where we have split (3.7a) and (3.7b) into free-stream terms, 

K@(u,u) := -{(D 2 u,D 2 ii) ,n + 2a 2 (Du,Du) ,n + ^(u,^^), (3.8a) 

K l b ° b ] (b,b) := -((D6,D6) ,fl + a 2 (b,b) , i2 ), (3.8b) 

contributions from the velocity profile U, 

K£3 (u, u) := -iaRe{{U D u, D u) , fl + a 2 (Uu, u) ,n - ((D D u)„,n), (3.9a) 

K^ 1 (6, 6) := -iaRm{Ub, 6) 0>fi , (3.9b) 

f roo _ ^xirf ^icc terms 

K[ s ]( U ,n):=-a 2 (S°( U )S^(n)*+S^( U )S°(«)*), (3.10) 
and contributions from the insulating boundary conditions 

K$(M) := -a(S?(6)S?(6)* + S°(&)S°(&)*). (3.11) 
Moreover, K„b : Vf, x 1— > C and K bu : V u x V b 1— > C are maps defined by 
K u6 (6, fi) = iaRe((B x Db,D fi) ,n + a 2 (B x b, u) ,« - ((D S x )6, D u) , fl ) 

- i?eB z ((D 6, D 2 fi)o, fl + a 2 (6, Du) ,n) - aRe S°(b)(iaB x (0) S° 2 (u) + B z S£(u))*, (3.12a) 
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K bu (u,b) = Rm{ia{B x u,b) ,n + B z (Du,b)o,n). (3.12b) 

For the parameterization (2.11) of the magnetic field we have 

K ub (b,u) = K l °l(b,u) + K l u B b \b,u) + K [ * b \b,u), Kbu(u,b) = K l b %,b) + K l b B J (u,b), (3.13) 

where, employing a similar notation as above, 

KJ2(6, u) := iaReA-W b, D u) 0i „ + a 2 (b, u) 0i „) - ReA~ 1 {{D b, D 2 u) ,n + a 2 {b, D u) ,n), (3.14a) 

Kfl(u,b) := iaRmA-\u,b) ,Q + RmA-^Du^n (3.14b) 

are free-stream terms, 

K !ri(M) := iaReH z Pm 1/2 ((BDb,Du) ,n + a 2 (Bb,ti) ,n - ((DS)6,Du) ,n), (3.15a) 

Kl?(M) := iaRmH z Pm 1 l 2 {Bu~b)vM (3.15b) 

are the contributions from the induced magnetic field B, and 

K S(M) := -aReSliWiaiA- 1 + A~ x RmB(0))S° 2 (u) + A^S^u))* (3.16) 

are free-surface terms. The maps K ua : C x 14 h> C, K ba : C x V b i— > C, and K aM : 14 x C ^ C, where 

K ua (a, fi) := -a 2 f — \ — + — ^ 2iaD £/(0)) aS^(u)* + ia(D 2 17(0) + H 2 D B(0))aS^(«)*, (3.17a) 



Pg l Re Oh 2 Re 

K ba {a~b) := iaA- 1 RmB B(0)aS^(b)\ (3.17b) 
K ou (u,a) :=S°(u)5*, (3.17c) 

represent the coupling of the velocity and magnetic field to the free-surface amplitude. Finally, Eq. (3.6b) 
contains the forms M uu : V u x V u i— ► C, : 14 x 14 i— > C, and M aa :CxC^C, where 

M uu (u,u) := Re((Du,Du) Q M + a 2 (B u,B u) a j 2 ), (3.18a) 



M bb (b,b) 
M aa (a,a) 



Rm(b,b) ,n, (3.18b) 
aa*. (3.18c) 

We are now ready to state the weak formulation of film MHD problems: 

Definition 1 (Film MHD problem). Let n = Q u V u = H 2 {f2), V b = H 1 ^), and V = V u x V b x C. Then, 
find (7,»)eCxl/\ {0}, such that for all v e V Eq. (3.5), with K and M given by (3.6), is satisfied. 

In a similar manner, one can construct weak formulations of the form (3.5) for channel MHD problems, as 
well as film and channel problems in the inductionless limit. In what follows, we will always use V to denote 
the full (direct sum) solution space. Also, we shall employ throughout the notation K and M for the stiffness 
and mass forms, and K UUl M uu , etc. for their constituent sub-maps. It is understood that the maps act on 
pairs of elements from the appropriate vector space, and their definitions are restricted to the problem type 
under consideration. In channel MHD problems, we select the solution spaces V u = Hq{Q), V b = i? 1 (J7), 
and V = V u x 14, where now fl = Q c . The stiffness and mass forms in (3.5) read 

K(v, v) = K uu {u, u) + K ub (b, u) + K bu (u, b) + K bb (b, b), (3.19a) 
M(v,v) = M uu (u,u) + M bb {b,b), (3.19b) 

where K bu , K bb , M uu , and M bb are defined as the corresponding maps for the film problem, i.e. (3.12b), 
(3.7b), (3.18a), and (3.18b). However, K uu and K ub are now given by 

K uu (u,u) = K{?](«,«) + («,«), K ub (b,u) = K [ °l(b,u) + K^ 1 (&,«), (3.20) 

where K„i, k[2, K^J, and K^j' have the same form as (3.8a), (3.9a), (3.14a) and (3.15a), respectively. 
The absence of the boundary terms in (3.20) is due to the essential imposition of the velocity boundary 
conditions (2.3). With these specifications, the variational formulations of channel MHD problems is as 
follows: 
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Definition 2 (Channel MHD problem). Set Q = Q c , and V u = H%{Q), V b = H 1 ^), V = V u x V b . Then, 
find (7, v) G C x V \ {0}, such that for all v G V Eq. (3.5), with K and M given by (3.19), holds. 

Film problems in the inductionless limit are governed by (2.2) subject to the boundary conditions (2.5), 
(2.7), (2.8) and (2.10). Like in full MHD problems we set V u = H 2 (Q), but now V = V u x C. The stiffness 
and mass forms then become 

K(v,v) = K uu (u,u) + K ua {a,u) + K au (u,a) + K aa (a,a), M(v,v) = M uu (u,u) + M aa (a,a), (3.21) 

where 

K uu (u,u) = K{?](u,«) + K{S(«,«) + ^l(u,u) + KW(u,«). (3.22) 
Here the form K^l : V u x V u 1— » C, defined by 

Kg4(u,ti) := -a 2 Hj(u,u) ,f2 + iair x H z ((Du,ti)o,n - (u,D«) 0l n) - H 2 (Du,Bu) ^, (3.23) 

represents the contributions from Lorcntz forces, and Kui, kL«', and Kui are given by (3.8a), (3.9a) 
and (3.10). Moreover, 

K ua (a,u) ~-o? ( + - 2iaDU(0)) aS° 2 (n)* - iaD 2 U(0)aSl(u)* (3.24) 

\Fg Re On Re J 

is the analog of (3.17a) in the inductionless limit. The maps K au , K aa , M uu , and M aa are defined in (3.17c), 
(3.7c), (3.18a), and (3.18c), respectively, and therefore we can now state the weak formulation of inductionless 
film problems: 

Definition 3 (Inductionless film problem). Let J? = J?f, V u = H 2 (f2), and V = V u x C. Find (-f,v) e 
CxV\ {0}, such that (3.5), with K and M given by (3.21), is satisfied for all v e V. 

The trial and test space for inductionless channel problems is simply V = V u = H 2 {f2). Moreover, the 
stiffness and mass forms reduce to 

K(v,v) = K uu (u,u), M(v,v) = M uu (u,u). (3.25) 

where 

K uu = K®(u,u) + K^(u,«) + KW(«,«), (3.26) 

and, as usual, Kuu, K^J and K^] are given by (3.8a), (3.9a) and (3.23), and M uu by (3.18a). Inductionless 
channel problems then have the following weak formulation: 

Definition 4 (Inductionless channel problem). Let V = H 2 (Q), where Q — Q c . Let also K and M be the 
stiffness and mass forms given by (3.25). Then, find v e V such that the relation (3.5) is satisfied for all v 
in V . 

We note that the weak formulation of non-MHD problems, in both film and channel geometries, follows 
by setting the Hartmann numbers Defs. 3 and 4 equal to zero. 

4. Galerkin Discretization 

The Galerkin discretization of the variational problems formulated in §3, collectively represented by equa- 
tions of the form (3.5), involves replacing the spaces V u and, where applicable, V b by finite-dimensional 
spaces V^" c V u and V b Nb C Vb, respectively of dimension N u and N b . Denoting the set of polynomials of 
degree p on Q by V p (f2), we define 

V?" := V u n T p "(f2), V b N " := V b n V Pb {(2), (4.1) 

where it is understood that fl stands for J?f (i? c ) when the problem under consideration is of film (channel) 
type. The subspaces V^ u and V b Nb provide a dense coverage of V u and V b in the limit N u , N b — > 00. 
Introducing the multi- index N, where N = (N u , N b ) for MHD problems, and N = N u for their inductionless 
counterparts, finite-dimensional spaces V N G V, where dimV N =: N can be constructed by substituting 
V^ u for V u , and V b Nb for V b in the definitions for V. Then, the Galerkin discretization of the variational 
problems (3.5) can be stated as follows: Find (7, v) € C x V N \ {0} such that for all v G V N the relation 

K{v,v)=')M{v,v) (4.2) 
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is satisfied. Here K : V N x V N ^ C is an approximation of K (the details of which will be specified in 
§4.2), oftentimes introduced to perform numerically the quadratures associated with the velocity and/or 
magnetic field profiles, U and B. However, in a number of cases, including the Poiseuille and Hartmann 
profiles considered below, the quadratures can be performed exactly and K = K for all elements of V-^. 
Given a basis {ipi}fLi of V-^, Eq. (4.2) is equivalent to the matrix generalized eigenproblem 

Kv = jMv, (4.3) 

where the stiffness and mass matrices, K g £ NxN and M g C NxN , have elements 

K mn = K(ip n , ip m ), M mn = M(ip n , ip m ), (4.4) 

and v — [v\, . . . , vn) t € C N is a column vector of the components of v in the {4>i}iLi basis. 

The matrix eigenproblem (4.3) can be solved using e.g. the QZ algorithm [47,48], or implicitly restarted 
Arnoldi methods [49]. Its numerical properties, such as roundoff sensitivity and memory requirements, 
depend strongly on the choice of basis for V^. Following the approach of KMS, in the sequel we use basis 
functions that are linear combinations of Legendre polynomials, constructed according to the smoothness of 
the underlying infinite-dimensional solution space (i.e. its Sobolev order), as well as the boundary conditions. 
In these bases, the matrices K and M are well conditioned at large spectral orders. Moreover, if the velocity 
and magnetic-field profiles are polynomial, they are banded and sparse, and stable, closed-form expressions 
exist for their nonzero elements. In problems with Hartmann steady-state profiles the stiffness matrix ceases 
to be sparse, and the contributions from U and B must be computed using numerical quadrature. Yet, K 
remains well-conditioned even at high spectral orders (see §5.2.2 below). 



4.1. Choice of Basis 



In order to construct our bases for and V h Nb we first introduce the reference interval fl := (—1,1) 
and the linear map Q : fl i—> fl = (z\, z 2 ), where 

Q(0 = zo+jti, z :={ Zl +z 2 )/2, j := h/2 := (z 2 - zf)/2. (4.5) 

In film problems we have z 2 — 0, z\ = —1, zq = 1/2, and h = 1, whereas in channel problems Q becomes the 
identity map (z 2 = 1, Z\ — —1, z = 0, and h = 2). Uniformly continuous functions on fl can be transported 
to fl via the pullback map Q* : C°(fl) ^ C a (fl), where (<3*/)(0 = f(Q(0) for an Y / e C°(Q). The 
pushforward map Q» : C°(fl) 1— » C°(fl), where (Q*f)(z) = f(Q^ 1 (z)) and / g C°(fl), carries out the 
reverse operation. Moreover, a straightforward application of the chain rule leads to the relations 

D dl (Q*/i)(Q(±l)) = r dl D d 7i(±l), (4.6a) 
((D d * 3 )D dl g»A,D d2 Qj 2 ) Q ,n = j 1 ^-^^ 9 (Q*g)t) dl A, f 2 ) a ^ (4-6b) 

where D is the derivative operator on fl, and and g are sufficiently smooth functions, respectively on fl 
and fl. 

Let L n , where n — 0, 1,2, . . ., denote the n-th Legendre polynomial defined on fl and normalized such 
that L n (l) = 1 (see e.g. [50] for various properties of the Legendre polynomials). The Legendre polynomials 
obey the orthogonality relation 

2J mn /(2n+l), (4.7a) 
where S mn is the Kronecker delta. In addition, the inner-product relations 



(wiL n ,L m ) Qi fi (n + l)5 m , n+ i n6. 



2 (2n+l)(2n + 3) (2n - l)(2n + 1) ' 

(w 2 L n ,L m ) Q fl (n + l)(n + 2)5 m ,n+2 . (n - l)(n + 1) + n(n + 2) n(n - l)S m , n - 2 



(4.7b) 



(2n+ l)(2n + 3)(2n + 5) (2n - l)(2n + l)(2n + 3) (2n - 3)(2n - l)(2n + 1) 

(4.7c) 
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hold, where Wk denote the weight functions Wfc(£) = £ fc - The values of the Legendre polynomials and their 
first derivatives at the domain boundaries are given by 

L„(-l) = (-1)", DL n (l)=n(n + l)/2, DL„(-1) = (-l)" +1 n(n + l)/2. (4.8) 

Moreover, the property 

(2n+l)L„ = D J L„-DL„_ 1 (4.9) 

is useful for evaluating integrals of the L n . 

We introduce the following linear combinations of Legendre polynomials, which will be used as bases of 
the vector spaces H${fl) n V P {Q) (for k G {0,1,2}), H 1 (Q) n V p (Q) , and H^(f2) nV p (f2): 
Proposition 1. The polynomials 



^•/,"* M > ^y ' (4 - iob » 



nV; 7_i ™ +lV7 ^2(2^ + 3) V 2n + 5 2n+l / ' v ; 

where 1 < n < N, span the spaces L 2 {fl)C\V N - 1 {fl), H^(f2)nT N+1 {f2), and H 2 { f2)r\T N+3 { f2), respectively. 
In addition, they satisfy the orthogonality relations 

(AM) 0>a = (DAW,DAW) 0/} = (D 2 aM,D 2 aSV, } = S mn . (4.11) 

Proposition 2. Let 

'(1-0/2, n=l, 

^(0:= (1 + 0/2, n = 2, (4.12) 
X-2(0. «>3. 

T/ien, {^n}^ =1 is a &asis of H 1 ^) n7 ,Ar_1 (/)). T/ie values of the basis functions at the domain boundaries 
are 

M-1)=M2(1) = 1, Mi(l) =M 2 (-1) = 0, (4.13a) 

Mn(±l) = 0, n > 3 (4.13b) 

Furthermore, the inner-product relations 

(D/ii,D/ii) 0iA = (DM2,LV 2 ) 0if } = -(D Ml ,D/x 2 ) 01 ^ - -1/2, (4.14a) 

(D/i n ,D// m ) ^ = 0, n e {1,2} and m > 3, (4.14b) 

(D/x n , DjU m ) 0j £ = (5 m „, n > 3 and m > 3 (4.14c) 

Proposition 3. TTie polynomials v n {Q, where 



"n(0 := < 



-(l + 0^-2)/4, n=l, 



(l + 2 (C-l)/4, n = 2, (4.15) 

IA[ 2 1 2 (C), 3<n<7V, 

span </ie space Hf(fi) C\V N+1 {f2). They have the properties 

^i(l)=LV 2 (l) = l, D^(1) = ^(1) = 0, (4.16a) 

i/„(-l) = LV„(-1) = 0, ne{l,2}, (4.16b) 

u n (±l) = Di/„(±1) = 0, n > 3, (4.16c) 
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and also satisfy the orthogonality relations 

(D 2 z. l7 D 2 ^) , } = -(B^,i)\ 2 ) 0/2 = 3/2, (D% 2 ,D% 2 ) f} = 2, (4.17a) 

^2 -^2 

(D i/„,D is m ) ft = 0, n e {1,2} and m > 3, (4.17b) 

- 2 - 2 

(D ^„,D v m ) j2 = dmn, n>3 and m > 3. (4.17c) 

One can check that within each of the sets of polynomials defined in Propositions 1-3 the elements are 
linearly independent and, as follows from (4.8), satisfy the appropriate boundary conditions. In particular, 
the polynomials A„ have the properties 

D J AW(±1) = 0, (4.18) 

where < j < k — 1 . In the context of FEMs they are referred to as internal shape functions of order k . On 
the other hand, /Ui and /x 2 , and v\ and v 2 are called nodal shape functions because they satisfy all but one 
of the conditions (4.18), respectively for k = 1 and k = 2. Separating the basis functions into internal and 
nodal ones facilitates the application of the natural boundary conditions at the free surface. For example, 
the forms (3.17) contribute only one nonzero matrix element, while (3.10), (3.11), and (3.16) contribute two. 

\k] . k 

Remark 5. The A„ polynomials embody H -regularity in the sense that they are, by construction, k-th 
antiderivatives of L 2 -orthonormal polynomials. As a result, the principal forms of the continuous spaces (i.e. 

* k * k ~ 

(D vi, D v 2 ) q for vi, V2 <G Hq(Q)) are, in accordance with (4.11), represented by identity matrices, and do 
not exhibit an element-growth problem with p. By virtue of (4.14) and (4.17), the corresponding matrices 
in the {/x„} and {f n } bases are, in each case, the direct sum of a 2 x 2 matrix and an identity matrix, and 
therefore are also well behaved. 

We obtain our basis polynomials for the discrete spaces and V b Nb (4.1) by transporting X„\ n n , and 
v n from the reference interval ft to the problem domain J? by means of the pushforward map Q*. Introducing 

^ _ [Q*v n film problems, (4 19) 

]<2*A^ channel problems, 

and Xn '■— Q*Hn, it follows from (4.1), in conjunction with Defs. 1-4, that V™ u = span{<3j n }^ 1 and 
V b Nb — span{xn}^=i- Then, our bases {ipn}n=i are constructed as follows: 

Definition 6 (Bases of the discrete solution spaces V^). In film and channel MHD problems we respectively 
set 

!(<£„, 0,0), l<n<N u , , 
(0,x»,0), N u + l<n<N u + N b , ^ := ' N ^ „ (4.20a) 

(0,0,1), n + N u + N b + l, l<°'*»>' ^ + l<n<N u + N b . 

Moreover, our basis vectors for inductionless film problems arc 

^ n :=(^' 0) < 1 (4.20b) 
\(0,1), n = N u + l, K ' 

whereas for inductionless channel problems we simply have ip n :— <p n , where 1 < n < N u . Thus, for all 
v G V N one can write v = X)^=ife]"V'n, where 

' (u T ,k T , a) , film MHD problems, 

(m t ,6 t ), channel MHD problems, (4 21) 

(u T ,a), inductionless film problems, 

m t , inductionless channel problems, 



with u e C N * andfteC^. 

We note here that the procedure of constructing finite-dimensional solution spaces by transporting polyno- 
mial functions from the reference element to the problem domain is extensively applied in ft,p-finitc-elcmcnt 
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Table 1 



Properties of the discrete spaces V u u and V b 

















Problem 


Definition 


Basis 


p 


Definition 


Basis 


P 


Film 
Channel 


Hl{Q t )nV p (fl{) 


{Q*Vn)n=l 


N u + 1 
Nu + 3 


^(flfJnFifl,) 


{<9*Mn}n=l 


N b - 1 
N b ~ 1 



methods (FEMs), with the difference that f2 is mapped to the mesh elements rather than the full domain 
J? (e.g. [27]). Our method can thus be viewed as a single-element hp-FEM, with h = 1 (film problems) or 
h = 2 (channel problems) . One of the benefits of working with afhne families of finite elements is that the 
action of sesquilinear forms on basis- function pairs only needs to be computed on (2, as the corresponding 
values on the mesh elements follow by scalings of the form (4.6). Even though this type of computational 
gain is not relevant to our single mesh-element scheme, working with £1 leads to a more unified treatment of 
channel and film problems, and also allows for extensions of the method to problems defined over multiple 
domains (e.g. vertically-stacked layers of fluids). The salient properties of the discrete solution spaces and 
their bases are displayed in Table 1. 

Remark 7. In channel MHD problems it is also possible to apply the procedure described by Shen [25,26] to 
construct linear combinations of Legendre polynomials that satisfy strongly (essentially) the Robin boundary 
conditions (2.4) for the magnetic field. That approach would lead to well-conditioned and (for polynomial U 
and B) sparse stiffness and mass matrices as well. However, since the boundary- value maps Sj(b) = D(b(zi)) 
cannot be defined for all elements of H 1 , the trial and test space for b would have to be an H 2 (D) subspace, 
such as Hl(Q) := {b e H 2 {Q); S\{b) - aSf(b) = S$(b) + aS$(b) = 0}. In our treatment of channel MHD 
problems, we opted to consider that b is an element of H l {Q) and enforce the boundary conditions weakly 
in the interests of commonality with our film-problem formulation. 

4.2. Structure of the Discrete Problems 

We are now ready to write down explicit expressions for the stiffness and mass matrices in (4.3). With 
the choice of basis functions in Def. 6, the free-stream contributions can be evaluated in closed form by 
means of the properties of the Legendre polynomials. This is also the case for the {/-dependent forms in 
problems with the Poiseuille velocity profile, since (4.7b) and (4.7c) can be used to evaluate terms that are, 
respectively, linear and quadratic in the reference coordinate £. On the other hand, the exponential terms 
in Hartmann profiles (2.12) preclude the derivation of closed-form expressions for the integrals, and one has 
to resort to numerical quadrature instead. Here we pursue two alternative approaches, either of which can 
be used to obtain highly accurate solutions of our stability problems. 

The first approach is based on specialized Gauss quadrature rules, by means of which the exponentially 
weighted sesquilinear forms are computed exactly (modulo roundoff error) . Numerical methods for orthogo- 
nal polynomials with exponential weight function over a finite interval, and the associated Gauss quadrature 
knots and weights, have been developed by Mach [34] . As with many classes of orthogonal polynomials, the 
challenge is to compute the coefficients of the three-term recurrence relation in a manner that is stable with 
the polynomial degree p. In the context of a study on optical scattering (a problem of seemingly little rele- 
vance to spectral methods), Mach presents an iterative algorithm that yields the required coefficients and, 
importantly, is stable at large p. By computing the eigenvalues and eigenvectors of the resulting Jacobian 
matrix (e.g. [51]), it is therefore possible to obtain quadrature knots and weights suitable for the evaluation 
of polynomial inner products weighted by exp(±iJ z z). 

We also propose an alternative approach, which, following the widely used practice in spectral methods 
( [35] and references therein) , involves replacing the weighted sesquilinear forms by approximate ones derived 
from numerical quadrature rules (in the present case, LGL quadrature). Banerjee and Osborn [37] have 
shown that in elliptical eigenvalue problems the incurred integration error docs not affect the exponential 
p-convergence of the discrete solution, provided that the eigenfunction being approximated is smooth, and 
the quadrature method is exact for polynomial integrands of degree 2p — 1. To our knowledge, however, 
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no corresponding theorem is available in the literature for the OS eigenvalue problems studied here, and, 
although surely an interesting direction for future research, an investigation along those lines lies beyond 
the scope of our work. Instead, in §5.2.2 below we contend ourselves with a series of comparisons with the 
exact-quadrature method supporting the adequacy of the 2p— 1 precision level in our schemes for free-surface 
MHD as well. 



4.2.1. Free-Stream Matrices 

For the matrix representations of the U and B-indepcndent forms it is convenient to introduce the square 
matrices T^t ld2 \ T^ 1 ^, anc [ yMi^ ^Qgg s j ze j s e q U£ j to the number of basis polynomials (i.e. N u 

and/or Nf,). Using, as in (4.7), Wk to denote the power-law weight functions Wk(Q — £ fc , we set 



-,[kd 1 d 2 ] 

ITT 

1 1 n 



:= KD d2 AW,D dl AM) 0/} , 



-,[kdid 2 
H? 



(w fc D 2 ^„,D 1 v m ) Q ft, 



-,[kdid 2 ] 
H 1 



(w k D 



(4.22a) 
(4.22b) 



For our purposes it suffices to restrict attention to the cases where all of k, r, di, and d 2 arc non-negative 
integers smaller than three. Then, closed-form expressions for (4.22), which we quote in Appendices A.l and 
A. 2, can be evaluated with the help of the orthogonality relations (4.11), (4.14), and (4.17). We note that 



several of the calculations can be performed in a hierarchical manner. Specifically, the property DA„ 
(see Proposition 1) carries over to the corresponding matrices, where the relation 



A 



n+l 



j-,[kdid 2 ] 






mn 



-,[fe,di-l,d 2 -l] 

h:- 1 



m+l,n+l 



(4.23) 



applies for r, d\, d 2 > 1. Moreover, by construction of the \i n and v n polynomials (Propositions 2 and 3), 
we have 



- 1 [kd 1 d 2 ] 
H 1 



-,[kdid 2 ] 



-,[kdid 2 ] 
H? 



-,[kd 1 d 2 ] 
H 2 



(4.24) 



where m, n > 3. That is, every N x N matrix j , ^ ld2 l contains a T^f ld ^ submatrix of size (N — 2) x (N — 2), 



and similarly a T^ dld2 ^ submatrix of size (N — 2) x (N — 2) is contained in every N x N matrix T^™ 

Remark 8. It follows from (4.7) that the matrices T^ ld2 ' are banded and sparse (see Table A.l). Moreover, 

their bands arc not fully populated, as every other diagonal consists of zeros. The bandwidth of T^t ld 
equal to 2r + k — d\ — d 2 . 

Remark 9. Let m and n respectively denote the row and column indices of T^ dld 

elements with m > 2 and n < 2, or m < 2 and n > 2, are the results of (weighted) inner products between 
nodal shape functions, respectively fii, [i 2 and v 2 , and the internal shape functions ^3, /U4, . . .and 1/3, 
V4, .... It can be checked by explicit calculation (see Appendix A. 2) that the spectral leakage between the 
nodal and internal shape functions is small. Specifically, the quantities 

l H 2 :=maxf M^ 1 ! ^ 0; n e {1, 2}) , l H , := max f \ T [ * dld2] ] ^0;ne{l,2}) (4.25) 

1 m 11 -"1 imn J m 11 n J mn J 

are found to obey the relation 



-,[kdid 2 ] 



is 



and T { ^ d2 \ Then, 



l H i = 4 + k — d\ — d 2 . 



(4.26) 



Note that defining l H 2 and l H i as maxima over the matrix columns n leads to the same expression as (4.26). 



In order to compute the matrix representations K^\, K^, and K [ ^u of the free-stream forms Kul (3.8a) 



(3.8b), and k£] (3.23), we employ the collective notation 



'[L] 



,[0] 



rj-i\kd\d 2 



' rp[kdid 2 ] 
rp[kd\d 2 \ 



film problems, 
channel problems, 



(4.27) 
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where T [kdld2] e R N »* N », and also write T l h k b dld2] := T^ ld2] e M^ xAr " in both channel and film problems. 
As above, we denote matrix rows and columns respectively by m and n. Then, making use of (4.6b), we 
obtain 



K [0] ._ 

-""66 ■ 



K 



^66 (Xm Xm) 



-J I J 



(> 



-2^(011] , ^2^(000] \ 
± 66 + " - 1 66 J 



and 



(4.28a) 
(4.28b) 

(4.29) 

(4.30a) 
(4.30b) 

The maps (3.14), coupling the velocity and magnetic fields, can be treated by introducing T^ kdld ^ € 

xAr ", where 



As for the mass forms M uu (3.18a) and M b& (3.18b), these are represented by the matrices 

M uu := [M uu (cp n ,cp m )} = Re j (j~ 2 T^ + a 2 T^ , 



M bb := [M 66 (x„,Xm)] = RmjT i 



[000] 
66 • 



xN " and Tg?^ ] G 



-,[fcditi2] 



:= (w fe D d2 A[ 2l ,D a Vm) 0/ }, 



rfi 



and also T 



[fcditZ 2 
bu 



mn 
xA S With 



-,[kdid2] 

mn 2 



:= (w k t) d2 v n ,T) dl Mm) ,fi) 



™[fedi<i2] 
- 1 6« 



-,[fcdid 2 ] 

' 

-,[fcdid 2 ] 



film problems, 
channel problems. 



Then, the matrices associated with k{^| and (3.14) are 



K bu - 



K ul{Xn,<j>m) 
K^J {4>n, Xm) 



iaReA, 



-,[011] 



iaRmA-^+RmA-^, 



- 2 tL7 ] ) 

[001] 



-,[021] 

ub 



+ a 2 T 



[010] \ 
«6 y 



(4.31) 

(4.32) 

(4.33a) 
(4.33b) 



where T [ ^ ] := (t£*" ,i] ) T 

4.2.2. {/ arte? B-Dependent Matrices 

We now examine the matrix representations of the forms k[S and K^' (3.9), and the maps k\^J and K^' 
(3.15), all of which involve inner products of Legendre polynomials with non-trivial weight functions. 

Problems with the Poiseuille profile (2.14) can be treated using the matrices T [kdld2] and T [ bb dld2] estab- 
lished in §4.2.1. First, we compute the action of the pullback map Q* (defined below (4.5)) on U, 

(Q*U)(0 = 1 - (Q(£)) 2 =: U + Utf + U 2 e, (4.34) 

where = 1 — z 2 , U\ = —ztyh, and U 2 = — h 2 /A. Specifically, in film problems (z\ = —1, z 2 = 0) we have 
Uo = 3/4, Ui = —1/2, and U 2 = —1/4, whereas in channel problems (z\ = —1, z 2 = 1) the trivial result 
Uo = 1, Ui = 0, and U 2 = —1 applies. We then set 



K [u] ._ 

-^66 ■ 



+ f/l [j- 
10] N 



---iaRej{Uo(3- 2 T^+a 2 T^ ] 
+U 2 {.r 2 T^+a 2 T^-2j- 2 T^ 



■ a 2 T [100] 



-2 T [010] 



(4.35a) 
(4.35b) 



where K [ %} g C JV « xJV " and K l b U b ] e C N " xN » respectively represent K^ 1 and K^ 1 . 
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Turning to problems with Hartmann profiles (2.12), it is convenient to introduce the shorthand notation 
= H z h/2, Sh ( {£) ■— s'mh(H^), and Ch ( {£) ■— cosh(_ff^^), which leads to the relations 



(Q*U)(0 = U - U s (0 - # c (0> (Q*B)(0 = -Bo - Bit + B s (0 + B c (0, 
with U = cosh(iJ z )/A, B = smh(H z )z /(H z X), B\ = smh(H z )h/(2H z X), and 

smh(H z z )s Hli j _ cosh(H z z )c Hi 6 _ cosh(H z z )s Hi ^ _ smh(H z z )c H( 



U s = 



B x 



B r 



(4.36) 



(4.37) 



X X H Z X ' c H Z X 

r rrl r rrl 

We also use £ l G k G (— 1, 1) and p G [, where H > and k € {1,2,..., G}, to denote the quadrature knots 
and weights computed via Mach's algorithm [34], such that 

d£e^/(0 = E4"L/O (4-38) 
fe=i 

holds for any polynomial / € 7 ,2G_1 (/?). Following the procedure outlined in Appendix A. 4, Eq. (4.38) can 

^« xA S Sjf ld2] G K^ xAr o, and 5[^ ld21 G K w » xA S where 



/: 



be used to evaluate the matrices S^ ld ' 2 ^ G 



..[deZicfe] 
'66 



..[deZicfe] 
*bu 



((D l> s )D "A^-D 'A^ 1 )^, channel problems, 
((D C/ S )D V„,D 1 ^m) o r }, film problems, 

d 2 



((D C/ S )D ' 

{((D -B S )D 2 A| 2 ],D Vm) /7, channel problems, 
((D B S )D 2 j/„,D Vm) ^ nml problems. 



[ddid2] 
66 



xAr \ and C 



[dd\d2 
ub 



(4.39a) 
(4.39b) 
(4.39c) 

»N u xN b 



Similarly, one can compute the matrices C { ^ lM G R^^", C 
whose elements are given by expressions analogous to (4.39), but with U s and B s respectively replaced by 
U c and B c . Then, K^} and become 



k'2 = -taRej (ft, (rML" 1 + a'Tir 1 ) - (r 2 sE" + a 2 <° 01 - r 2 s£°i) 



(> 



4 C 



[on] 



v 2^[000] 



•r 2 < 01 



= —\aRmj (Uqj 



,[000] 



,[ooo] 



<° 01 ) • 



(4.40a) 
(4.40b) 



Also, the maps K^ 1 and K'f 1 (3.15) induce the matrices K [B] G C JV « xJV » and K [ b B J G C" 1 ^" given by 



K [B] ._ 



iaReH z Pm 1/2 j 
-Bi (j 



-Bo (./ ' • 



2^.[011] 

i/6 



2c? l000] ,-2 c [110] 



[000] 
ub ^ a ± ub 



, -oL- 



[000] \ 

■,[011] 

' «6 



i/^L«^j , ^2^(000] •-2 r ,[H0] 



wb 



2y.[ m ] 



[010] \\ 



= mllmlLPmWj (-B T™ B,T^ + S f^ + Cf^) 

T 



where := (s 



[dd 2 di 
bu 



ub )) 

[10( 

bu 



(4.41a) 
(4.41b) 



and C|l dlti2] := (c [ £ ld ' ] ) 

Remark 10. Due to the non-polynomial nature of the Hartmann profiles (2.12), the matrices in (4.40) 
and (4.41) are fully populated, and no simple closed form expressions exist for their evaluation (cf. (4.35)). 
However, by virtue of (4.38) no quadrature errors are made in the computation of their elements. 

The expressions presented thus far are restricted to the specific cases of the Poiseuille and Hartmann 
profiles. Oftentimes, however, one is faced with the task of studying the stability properties of arbitrary 
steady-state profiles, and, although in principle possible, deriving each time specialized quadrature schemes 
would be a laborious task. An alternative approach is to replace the weighted forms and maps by approximate 
ones defined on the discrete solution spaces by means of the following procedure: 
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Let (o,k € [z\, z 2 ], where k = 0, 1, . . . , G + 1, Co = z i, and Cg+i = z 2, be the abssicas of LGL quadrature 
with G interior points on the interval D — [z\, z 2 \, and let Q Gyk be the corresponding weights (this type of 
quadrature is exact for polynomial integrands of degree up to 2G+1 [51]). Also, consider inner products of the 
form (Wfi,f2)o,o, where W stands for either U or B, or their derivatives, and /1 and f 2 are polynomials 
of degree p\ and p 2 , respectively. For all such inner products appearing in K^t, KjV', Kj?, and K 



set 



G > {pi + p 2 )/2 - 1 and make the substitution {Wf 1 ,f 2 )o,n ^ J2k=o QG,kW{^ G , k )h{C,G,k)fc(fe,k)*- The 
resulting forms and maps, respectively denoted by K 



V., 



N„ 



x V 



K 



Kit 1 : V b N » x V"- 



C, and k[ B J ■ 



bb 



V h N " x V"* 



: V™ u x Vl 



C, are 



KM(u,u) :=-iaRe £ g G , k (U(( G , k )(D u(C G , fe )D u(( G , k )* + a 2 u(C G , k )u*(( G , k )) 

k=0 

-(DC/(C G>fe ))n(C G , fe )Dn(C G , fe )*), 

G b +1 

Kf b \b,b) := -iaRm ]T Q G , k U{C G , k )b{C G , k )b(C G , k )\ 



(4.42a) 
(4.42b) 



fe=0 



and 



where 



G„ b +i 

K [ *\b,u) ^laReHzPm 1 ' 2 £ Q Gttb<k (B(( Gub , k )(D b(( Gub , k )D u(( G ^ k )* 

fe=0 

+a 2 b(( Gub , k )u(( Gub , k y) - (D S(C Gut ,fc))6(CG ut ,fc)D fi(CG ut ,fc)*) , (4.43a) 

G„„ + l 

K[ B u ] (u,b) :=\ a RmH z Pm 1 / 2 £ ^.^(Cg^/O^Cg^WCg,,,,,*)*, (4.43b) 

fc=0 

G M > p u - 1, G 6 > P6 - 1, G u() > (p u + p 6 )/2 - 1, (4.44) 

and, as usual, p u and p b are the polynomial degrees of the velocity and magnetic-field bases (see Table 1). 
The sesquilinear form K : V x V N 1— ► C introduced in (4.2) then follows by replacing the exact forms and 
maps in (3.6a) with the corresponding approximate ones defined in (4.42) and (4.43). 

Remark 11. Our choice of precision in (4.44) is motivated by Banerjee and Osborn's [37] result that in 
finite-element methods for elliptical eigenvalue problems it suffices to use quadrature schemes that are exact 
for polynomial integrands of degree 2p — 1, where p is the degree of the FEM basis. Here we do not pursue 
a formal proof of the adequacy of (4.44), but the numerical tests in §5.2.2 demonstrate that eigenvalues 
computed using the smallest quadrature precision consistent with it converge in a virtually identical manner 
with those obtained via the exact quadrature scheme. 

For the purpose of evaluating the matrices representing , k[^' , R[f , and Rj^' , which we again denote 

G b xN b 



,xJV„ 



bv K [u] K [U] 



K^ b \ and K b ^} , we introduce the differentiation matrices A$ G M GuXNu 



kn 



and A' b [d] G R c 

\l) d is n (( Gu . k ), 
lD d AL 21 (CG„,fe), 



,xJV, 



6 with elements 



kn 



film prob., 
channel prob. 

:= D fJ, n (CG b ,k), 



kn 



< d] 



kn 



|D v n {Ca ub , k ), 

:= D Hn{C,G ub ,k), 



film prob. 
channel prob. 



(4.45a) 



(4.45b) 



where Q G , k G [— 1, 1] for k G {0, 1, . . . , G+ 1} are LGL quadrature knots on the reference interval fl. We also 
make use of the G x G diagonal weight matrices q g , whose entries [f? G ]fcfc = Q G , k are equal to the quadrature 
weights associated with the knots ( Gjk (note that Q Gtk = Q~ 1 {( G , k ) and g Gjk = 2g Gyk /h), and the diagonal 



matrices U^q and Bg, where 



U 



kk 



:=D Q*(U)(Ca,k) and 



B 



[d] 



kk 



:= D Q*{B)(( Gtk ). We then obtain 



*E3 ■= [^(<PnM] =-iaR ej (r^YeGU§ u A^ + a 2 (A^QGM[^ ] 
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K 



[U] 



bb 



^bb (Xn,Xm) 



-iaRm 3{ ^fg Gb U^Af, 



(4.46a) 
(4.46b) 



and 



K 



IB] ._ 



ub 



^bu - 



^ [ ub\Xn,<t>m) 



^bu (4>n,Xn 



\ a ReH z Pm>/ 2 j{j- 2 {Af ] ) T QG u 



(-rub O 



laR m H z P m ^j{Af^g Gub B^AT 



[0] 
6 

(4.47a) 
(4.47b) 



Remark 12. In 64-bit arithmetic the numerators and denominators in (2.12) overflow at around H z = 
ln(2 1023 ) ~ 700 . This, in conjunction with the fact that neither U nor B have Taylor expansions about 
H z = oo valid for all z € [—1,1], renders the evaluation of the U and B matrices at large H z somewhat 
problematic. A practical workaround is to code the internal calculations for U and B in REAL*16 (128-bit) 
arithmetic, supported by a number of Fortran compilers (e.g. the Intel and NAG compilers), pushing the 
occurrence of the overflow to H z = ln(2 16 ' 383 ) ~ 11,000. Note that a similar issue arises with the exact- 
quadrature method, but in that case performing the internal computations with REAL*16 data types is not 
as straightforward (see Remark 16 in Appendix A. 4). 



4.2.3. Boundary Terms 

Boundary terms, namely (3.10), (3.11), (3.16), (3.17), and (3.24), are the outcome of the natural imposition 
of the stress and kinematic conditions at the free surface, and the insulating boundary conditions at the wall 
and the free surface. One of the benefits of working in the {fi n }n=i an d { v n}n=i bases, consisting of internal 
and nodal shape functions, is that each of the boundary terms contributes at most two, p-independent, 
nonzero matrix elements in the stiffness matrix K. In consequence, its sparsity and conditioning are not 
affected by the boundary conditions. Specifically, Kul, K^, and K^j are represented by the matrices A"|fJ e 



xA \ and K [ *l e C N » xNb , where, as follows from (4.13) and (4.16), 



K [S] 
K, 



" UU 

'[i] 



'bb 
K ub 



ni Xm, 



^ubiXn,<t>m) 



= ~j 1 a 2 (5 m iS n2 + 5 m2 S n i), 

= aRe (iaiA- 1 + A' 1 RmB(0))6 ml - j^A^S^) 8 n2 . 



In addition, the maps K ua , Kfc a , and K au respectively give rise to the column vectors K_ v 
K hn e C Nb , and the row vector K J„ e M. Nu , where 



e C 



(4.48a) 
(4.48b) 
(4.48c) 
" and 



K UO (0„,1) = a 2 



-1 



a 



Pg 2 Re 



Oh 2 Re 



[KbJn ~ Kba(Xn,l) = iaAj 1 RmDB(0)S n2 , 



+ 2ia D U(0) )S nl + j (D 2 17(0) + H 2 D B(0)) S n2 , (4.49a) 

(4.49b) 
(4.49c) 



In inductionless problems the column vector corresponding to (3.24) is given by (4.49a) with D B(0) formally 
set to zero. 



4.2.4. Constructing the Stiffness and Mass Matrices 

Eq. (4.4), according to which the stiffness and mass matrices are to be computed, has different instantia- 
tions, depending on the forms K and M of the variational problem at hand (Defs. 1-4), and the corresponding 
choice of basis functions (Def. 6). The matrices introduced in §4.2.1— §4.2.3 serve as building blocks, out of 
which K and M can be composed in a modular manner. A number of these matrix 'modules' are common 
among different types of problems (e.g. (4.28a) and M uu (4.30a) are present in all film and channel 

problems), which is convenient for implementation purposes. 
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In film MHD problems, K and M have N — N u + N b + 1 rows and columns, and are given by 



K 



K 



UU 

I K 



K 

[U] 

UU 




K 



[o] 
ub 



K 



K 

[S] 



ub 



M 



M uu 
M bb 
1 







(4.50) 



*U(0) J 



Here the submatrices with uu and bb indices are respectively dimensioned N u x N u and N b x N b . Among 
them, the [/-independent matrices, K^\, K^\, K b ° b , K b \, M uu , and Mbb, are given by Eqs. (4.28), (4.48a), 
(4.48b), and (4.30). Also, the submatrix K^} is to be evaluated using either of (4.35a), (4.40a), and (4.46a), 
depending on whether the velocity profile is Poiseuillc, Hartmann (treated by means of the exact-quadrature 
method), or LGL quadrature is employed. Similarly, K b ^ can be computed by means of either (4.35b), 
(4.40b), or (4.46b). The submatrices indexed by ub and bu have dimension N u x N b and N b x N u , respectively. 



K^ b \ and K^, follow from (4.33) and (4.48c), while for those that depend 

-IB] 



The ^-independent ones, K 

on the induced magnetic field, namely and K^ 1 , there exist options to use exact quadrature (4.41), or 
LGL quadrature (4.47). Finally, the column vectors K_ ua and K_ ba , respectively of size N u and N b , and the 
row vector K^ lt of size N u are given by (4.49). In inductionless film problems (Def. 3), the magnetic-field 
degrees of freedom are not present, and K and B are replaced by the (N u + 1) x (A^„ + 1) matrices 



(4.51) 





( + K w 

UU 1 UU 
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K = 


+K [u\ + K m 




, Af = diag(Af „„,!), 




{ Elu 


-ia?7(0) J 





where, aside from which is given by (4.29), and K_ ua (obtained from (4.49a) with DB(0) set to zero), 

the submatrices have the same definitions as in (4.50). 

In the interests of brevity, we do not write down explicit expressions for the stiffness and mass matrices in 
channel problems. We note, however, that they have the same structure as the corresponding film-problem 
matrices, but with the rows and columns representing the free-surface removed, and all boundary terms 
involving the velocity {K^\, and K^ b ) set to zero. 

Remark 13. The mass matrices in (4.50) and (4.51), as well as in the corresponding channel problems, 
are symmetric positive definite (SPD). Rewriting (4.3) in the form -jmKv = jkMv, where jk/jm = 7 
(the QZ algorithm [47] actually solves this version of the problem), the non-singularity of M guarantees 
that 7m ^ (i.e. 7 is finite). In fact, as can be checked from (3.6b), M is SPD for all choices of discrete 
bases. In tau methods, however, M can be singular. Dawkins et al. [24] have shown that in the Legendre 
tau discretization of a fourth-order eigenvalue problem (structurally similar to the OS equation) M has a 
non-trivial nullspace, and, as a result, the discrete problem contains an infinite eigenvalue. Moreover, the 
Chebyshev tau formulation of the same problem was found to contain spurious eigenvalues, even though 
in that case M is non-singular. Treating the Legendre and Chebyshev tau methods as members of the 
one-parameter family of Gegenbauer tau methods, the spurious eigenvalues in the Chebyshev case were 
interpreted as perturbations of the infinite eigenvalues in the Legendre one. Like KMS, we found no evidence 
of spurious eigenvalues in any of the schemes presented here, which, in light of the analysis by Dawkins et al., 
is likely due to the fact that the variational formulation described in §3 leads to non-singular mass matrices 
irrespective of the choice of basis. 

Remark 14. The sparsity of K and M in problems with polynomial steady-state profiles enables the 
efficient use of iterative solvers. A number of implementations (e.g. the ARPACK library [49], which is 
also available in Matlab) provide the option to specifically seek the eigenvalues with the largest real parts, 
which are oftentimes the ones of interest. In practice, however, we observed that these are particularly hard 
eigenvalues to compute, with the algorithm frequently failing to achieve convergence. Instead, we found that 
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a more feasible strategy is to search for eigenvalues with the smallest absolute value. Due to the predominance 
of highly damped modes in the spectrum (i.e. eigenvalues with large I7I but small Re(7)), the eigenvalue 
with the largest real part often happens to be among the smallest absolute value ones. This approach was 
used to compute the eigenvalues at p = O(10 3 ) in Fig. 11 below. 

5. Results and Discussion 

In this section we present a series of test calculations aiming to validate our Galerkin schemes, and 
illustrate the basic properties of our stability problems. First, in §5.1 we study the eigenvalue spectra of 
representative film and channel problems. Various aspects of numerical accuracy are examined in §5.2. In 
§5.3 we test the consistency of our schemes against energy conservation in free-surface MHD, and the time 
evolution of small-amplitude perturbations in nonlinear simulations. The critical-parameter calculations in 
§5.4 is our final topic. Aside from the nonlinear simulations in §5.3, all numerical work was carried out using 
a Matlab code, available upon request from the corresponding author. We remark that in order to facilitate 
comparison with relevant references in the literature, we oftentimes express our computed eigenvalues in 
terms of the complex phase velocity c = ry/a, rather than the complex growth rate 7. As stated in §2.1, in 
the former representation a mode is unstable if r := Im(c) > 0, while C :— Re(c) its phase velocity. 

5.1. Eigenvalue Spectra of Selected Problems 

5.1.1. Non-MHD Problems 

One of the most extensively studied problems in hydrodynamic stability is non-MHD channel flow with the 
Poiseuille velocity profile (e.g. [6,7,11,12,21,52]). In the high Reynolds number regime, the spectrum of the OS 
operator forms three branches on the complex plane, conventionally labeled A, P, and S [32]. The branches 
are shown in the numerical spectrum in Fig. 3, obtained at Re = 10 4 and a — 1 by solving the matrix 
eigenproblem (4.3) derived from Def. 4 (with H x = H z = 0). According to Table 1, u{z) is expanded in the 
basis, where for the present calculation the polynomial degree is set at p u = N u + 3 = 500. Due to 
the reflection symmetry of (2.2) and (2.3) with respect to z, the cigenfunctions fall into even (u(—z) = u{z)) 
and odd (u(— z) = —u{z)) symmetry classes. The S branch contains a countably infinite set of modes, whose 
phase velocity is asymptotically equal (at large and negative Im(c)) to the mean basic flow (U) = 2/3 (2.13). 
On the other hand, the A and P branches contain a finite set of modes, respectively with < C < (U) and 
(U) < C < 1. The P modes come into nearly degenerate even and odd pairs. As noted by Orszag [11], this 
near degeneracy is a genuine property of the spectrum, which does not disappear by increasing the spectral 
order. While all of the P modes are stable, the A branch contains a single unstable mode of even symmetry. 
This instability is of the critical-layer type [21]: At sufficiently high Reynolds numbers, and over a suitable 
range of wavenumbers, the energy transfer from the basic flow to the mode (the Reynolds stress (2.16a)) 
exceeds the viscous dissipation, and as a result its growth rate becomes positive. 

Table B.l lists in order of decreasing Im(c) the first 33 eigenvalues plotted in Fig. 3. This calculation 
has previously been performed by Kirchner (see Table VII in [6]) using the same Galerkin scheme as in the 
present work, so the two sets of eigenvalues should be in very close agreement. A comparison (see also the 
underlined digits in Table B.l) reveals that for modes at the top end of the spectrum the relative agreement 
is of order 10 -15 , i.e. close to machine precision. However, descending the spectrum, the number of decimal 
digits for which the calculations agree exhibits a diminishing trend, culminating to an O(10 -9 ) relative 
difference for Mode 33. This discrepancy is likely due to roundoff sensitivity in the computed eigenvalues 
close to the intersection point between the A, P and S branches, which is known to increase steeply with Re 
[33]. In our schemes, machine roundoff in double- precision (64 bit) arithmetic already leads to relative errors 
of order unity at Re <~ 5 x 10 4 (see §5.2.3 below). Therefore, the observed six-digit loss in the agreement 
between Kirchner's eigenvalues and ours is not unreasonable at Re = 10 4 , especially for modes like A10, 
which lies particularly close to the intersection point (Im(c) = 0.637 ~ 2/3). 

In film problems, again with the Poiseuille profile, the eigenproblem (4.3) is derived from Def. 3 (with 
H x = H z = 0), and, in accordance with Table 1, the velocity eigenfunction is expressed in terms of the v n 
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Fig. 3. Spectrum of a non-MHD channel problem with the Poiseuille velocity profile at Re = 10 4 , a = 1, and p u = 500, showing 
the A, P, and S branches, o and X markers respectively correspond to even and odd modes. The even mode marked in boldface 
is unstable. 

polynomials. Our nominal specification of the free-surface parameters (which will also be used in the MHD 
calculations below) is Oh — 3.14 x 10~ 4 and Pg = 1.1CP 4 , corresponding to a typical liquid-metal film of 
thickness 0.01 m at terrestrial gravitational fields [8]. Setting a = 1 and p u = N u + 1 = 500, we evaluate 
the spectra at Reynolds numbers Re = 10 4 and Re = 3 x 10 4 . The resulting eigenvalues are displayed 
in Fig. 4 and tabulated in Table B.2, which also lists the modal free-surface energy (2.15c). Like channel 
problems, the spectra exhibit the A, P, and S branches, and additionally contain two modes associated 
with the free surface, labeled U and F. Mode F is a 'fast' downstream-propagating surface wave, whose 
phase velocity is always greater than the basic velocity at the free surface (Re(c) > 1). It is unstable for 
Re > (5/8) 1 / 2 /P 5 [8,53], provided that a is smaller than some upper bound. This so-called soft instability is 
present in Fig. 4(b). Mode U is an upstream-propagating mode (Re(c) < 0), which is present in the spectrum 
at sufficiently small Reynolds numbers (e.g. Fig. 4(a)). For Re < 10 3 (and a = 1) its eigenfunction has the 
characteristic exponential-like profile of a surface wave. However, as Re grows its phase velocity increases, 
because the mode tends to be advected downstream by the basic flow. At the same time, its eigenfunction 
develops the characteristics of an internal (shear) wave, such as well defined boundary and critical layers. 
Eventually, the eigenvalue crosses the Re(c) = axis and merges with the A branch, taking over the role 
of the Ai mode in channel flow (for this reason, in Table B.2 Mode U is also labeled At). In particular, 
provided that the Reynolds number exceeds some critical value, it experiences an instability very similar 
to that in channel flow, oftentimes referred to as the hard instability. The spectra in the top and bottom 
panels of Fig. 4 respectively lie below and above the hard-instability threshold. As can be checked from 
Table B.2, only a relatively small number of modes carry appreciable free-surface energy. Apart from the F 
mode, and the upper A and P family modes, for which E a /E <~ 0.5, the remaining modes are internal, with 
E a /E<\{r\ 

5.1.2. Problems in the Inductionless Limit 

The simplest version of MHD is the inductionless approximation (2.2), whose weak formulation is stated 
in Defs. 3 and 4, respectively for film and channel problems. Compared to the non-MHD baseline scenario, 
the steady-state magnetic field, parameterized by the streamwise and flow-normal Hartmann numbers H x 
and H z , affects the stability of the flow both at the level of the basic state, as well as the perturbations. In 
the former case, the flow-normal component of the field leads to the establishment of the Hartmann velocity 
profile (2.12), which differs substantially from the Poiseuille one, even at moderate Hartmann numbers (see 
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Fig. 4. Eigenvalues of non- MHD film flow with the Poiseuille velocity profile at Oh = 3.14 X 10~ 4 , Pg = 1.10 X 10~ 4 , a = 1, and 
p u = 500. The Reynolds numbers are Re = 10 4 (a) and Re = 3 X 10 4 (b). In addition to the A, P, and S branches encountered 
in channel problems (Fig. 3), the spectra contain downstream-propagating surface waves with Re(c) > 1, labeled F. Also, an 
upstream-propagating wave (Re(c) < 0), labeled U, is present in the Re = 10 4 spectrum. At Re = 10 4 all of the modes have 
negative growth rates, but at Re = 3 X 10 4 Modes Ai and F (represented by boldface markers) are unstable (Im(c) > 0). 



Fig. 2). The departure from the parabolic profile affects the Reynolds stress, which is the main driver of 
critical-layer instabilities. The magnetic field also acts at the level of the perturbations by way of electrical 
currents induced by the perturbed fluid motion within the field. These induced currents set up Lorentz 
forces, which, in accordance with Lentz's law, always tend to dampen the flow. Moreover, they modify the 
velocity distribution of the perturbations, changing in turn the Reynolds stress and/or viscous dissipation. 
It is generally known, both on theoretical grounds [54,55], as well as from numerical calculations [16,19], 
that in channel problems the combined outcome of these effects is strongly stabilizing. In film problems, 
however, the existence of a resonance between the velocity and surface degrees of freedom may cause the F 
mode to deviate from that behavior [8]. 

We first consider film problems with flow- normal magnetic field (H x = 0). Fig. 5 displays the eigenvalues 
computed at H z = 14 and 100, with all other spectral and flow parameters equal to those in Fig. 4(b). 
Numerical results obtained using both exact and LGL quadrature for the computation of the stiffness matrix 
K (respectively (4.40a) and (4.42a)) are listed in Table B.3. The maximum relative difference between the 
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Fig. 5. Eigenvalues of inductionless film problems with the Hartmann velocity profile and flow-normal background magnetic 
field (H x = 0) at Re = 3 X 10 4 , Oh = 3.14 X 1CT 4 , Pg = 1.10 X 10~ 4 , a = 1, and p u = 500. The flow-normal Hartmann 
numbers are H z = 14 (a) and H z = 100 (b). 

two eigenvalue sets is of order 10~ n (the corresponding mode is P22 at H z = 14), becoming as small as 
4.5 x 10~ 16 for the P3 mode. We note that the agreement between the lower modes does not improve by 
increasing p u . Like in our previous comparison of the eigenvalues of plane Poiscuillc flow (Table B.l) with 
the corresponding calculations by Kirchner [6] , the numerical convergence of the lower modes appears to be 
over fewer significant digits than the least stable ones. Nonetheless, our results demonstrate that the LGL 
quadrature scheme is a very viable alternative to the exact one, especially in light of its flexibility to treat 
arbitrary analytic velocity profiles (see also §5.2.2 ahead). 

Comparing Fig. 5 to Fig. 4(b) illustrates the following basic aspects of the magnetic field's influence on the 
eigenvalues. First, as H z increases the A branch is seen to collapse. That is, the eigenvalues move towards 
the intersection point between the P and S branches, eventually experiencing what qualitatively appears as 
an inelastic collision with the S branch. In the process, Mode Ai (the hard mode) crosses the Im(c) = axis, 
i.e. is stabilized. The real part of the S family eigenvalues remains (asymptotically) equal to the average 
value of the velocity profile, and moves from 2/3 towards 1, in accordance with (2.13). At the same time, 
the P branch becomes progressively aligned with the S branch. For sufficiently small values of H z , including 
the H z = 14 example in Table B.3, the P modes are somewhat less stable than in the non-MHD case 
(cf. Table B.2), but never cross the Im(c) = axis. As for the originally unstable F mode, this also becomes 
stabilized once H z exceeds some critical value (the spectra in Fig. 5 and Table B.3 are evaluated past that 
threshold) . 

The behavior outlined above is encountered at moderate H z , and is mainly due to the formation of the 
Hartmann velocity profile. As discussed in Rcf . [8] , at sufficiently large Hartmann numbers Lorentz damping 
causes the decay rate |Re(7)| of the A, P, and S modes to increase quadratically with H z . In contrast, 

1/9 

once H z crosses a threshold scaling like Pg 1 , the F mode exhibits a change in behavior, with its decay 
rate switching over to a decreasing function of the Hartmann number. In the H z = 100 problem in Fig. 5 
and Table B.3, which lies close to that transition, the decay rate \r\ — 0.12765 of the F mode already is 
substantially smaller than that of the Lorentz-damped P and S modes (\r\ > 0.31013). A single A mode is 
present in the spectrum with comparable decay rate \r\ = 0.13099, but at larger Hartmann numbers (not 
shown here) it also becomes suppressed. Even though the F mode remains stable, eventually it becomes 
the only one with small decay rate. As a result, film and channel problems differ qualitatively in that the 
eigenmodes of the former cannot be damped by arbitrarily large amounts solely by applying a background 
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Fig. 6. Spectrum of inductionlcss channel flow with the Hartmann velocity profile and oblique background magnetic field 
(orientation angle 4> = 1° with respect to the streamwise direction) at Re = 10 4 , H z = 14, H x = 14/tan(</>) 802, a = 1 and 
Pu = 500 

magnetic field. 

The more general case with oblique external magnetic field (i.e. H x and H z both nonzero), is especially 
interesting in the context of channel problems because, as can be checked from (2.2), the reflection symmetry 
with respect to z is no longer present. As shown in Fig. 6, the near-degeneracy between the even and odd 
P-family modes is broken, and the resulting sub-branches assume a distinctive curved shape. In general, 
the streamwise Hartmann number required to cause a comparable change in the eigenvalues is significantly 
larger than the corresponding flow-normal one. It is for this reason than in Fig. 6 we consider a magnetic 
field oriented at only 1° with respect to the streamwise direction, but sufficiently strong so that H z is equal 
to the one used in Fig. 5(b). In film problems, where no nearly degenerate P mode pairs exist to begin with, 
the oblique field still causes the P branch to adopt a qualitatively similar curved shape. 

5.1.3. Film MHD Problems 

We now relax the inductionless approximation made in the preceding section and consider film MHD 
problems, defined variationally in Def. 1, and discretized using the {fi n }n=i an d { v n}n=i bases for the 
magnetic- field and velocity eigenfunctions (see Table 1). Throughout this section we work at polynomial 
degrees p u = Pb = 500. Moreover, we compute the U and independent terms in the stiffness matrix K using 
the exact quadrature scheme, i.e. (4.40) and (4.41), although accurate results can also be obtained by means 
of the LGL method (Eqs. (4.46) and (4.47)). 

Noting that the limit Pm — > 0, at which Eqs. (2.1) reduce to (2.2) (under the proviso that H x and H z are 
non-negligible), is a singular limit of the coupled OS and induction equations, one can deduce that certain 
MHD modes, which we refer to as magnetic modes, are disconnected from the inductionless spectra. These 
are to be distinguished from hydrodynamic modes, that are regular as Pm — > 0. Whenever Pm is of order 
unity, magnetic modes are expected to be present in the portion of the complex plane with Im(c) > — 1, 
irrespective of the background magnetic-field strength. In fact, they stand out particularly clearly in spectra 
evaluated at H x — H z = 0, such as the one depicted in Fig. 7 for a Pm = 1.2 problem. In this special case 
with zero background field, the maps K„6 (3.12a), Kb u (3.12b), and Kb a (3.17b) vanish, and the magnetic 
modes are independent of the hydrodynamic ones. The latter have zero magnetic-field eigenfunction and 
retain the same velocity eigenfunction and free-surface amplitude as in the non-MHD case, whereas for the 
former u and a are zero and b is non-vanishing. The magnetic modes form a three-branch structure as well, 
whose branches we label Am, Pm, and Sm. The magnetic S branch coincides with the hydrodynamic one, 
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Fig. 7. Eigenvalues of film MHD flow with the Poiseuille velocity profile and vanishing steady-state magnetic field (H x = H z = 0) 
at Re = 10 4 , Pm = 1.2, Oh = 3.14 X 1(T 4 , Pg = 1.10 X 1CT 4 , a = 1, and p u = Pb = 500. In addition to the hydrodynamic 
modes, marked with o, the spectrum contains magnetic modes, labeled by +, which form the Am, Pm, and Sm branches. 

and the Pm branch lies close, but does not coincide, with P. The Am branch forms a nearly straight line 
that interpolates between the hydrodynamic A modes. Numerical values for the complex phase velocities of 
the 25 least stable magnetic modes are tabulated in Table B.5. 

When the steady-state magnetic field is nonzero, Kb u , ^ub, and Kb a couple the hydrodynamic and mag- 
netic modes, typically resulting to the formation of multiple eigenvalue branches. This type of behavior is 
illustrated in Figs. 8 and 9 for film MHD problems at Pm = 1.2, respectively with flow-normal and oblique 
external magnetic field. Tables B.6 and B.7 list the corresponding complex phase velocities and energies. As 
shown in Fig. 8(a), instead of leading to the collapse of the A branch and alignment of the P and S branches 
observed in the inductionless limit (Fig. 5), the magnetic field causes the nearly coincident three-branch 
structures at H z = to split into two distinct ones, each of which is populated by both hydrodynamic 
and magnetic modes. Moreover, an unstable magnetic mode (Mi) is now present in the spectrum. This 
mode, which also arises in channel problems, signifies that at sufficiently high Pm the magnetic field can 
destabilize an originally stable flow. As H z increases above 14, the tails of the branches split again, result- 
ing to the intricate eigenvalue distribution observed at H z = 100, which, apart from Mode Mi, is nearly 
symmetrical about Re(c) = 1. The spectrum with oblique external magnetic field (Fig. 9) exhibits multiple 
branches as well, and additionally contains a second unstable magnetic mode (M2). In both examples with 
H z = 100, Mode M 2 stands out in that its kinetic energy is significantly smaller (E u /E = 0.0023948 and 
E u /E = 0.00091314 for H x = and H x = if 2 /tan(l°), respectively) than the E u /E = O(10" 1 ) values of 
the remaining modes. 

At the Pm < 10~ 4 regime of laboratory fluids, we have observed that film MHD spectra are well approx- 
imated by those in the inductionless limit, apart from the presence of (i) magnetic modes with Im(c) > — 1, 
and (ii) a diffusive interaction between the F and P modes, accompanied by an instability. These two kinds 
of discrepancy are shown in Fig. 10 and Table B.8, where Pm = 10~ 4 and inductionless spectra have been 
evaluated at Re = 10 6 , a = 0.01, {H X ,H Z ) = (0, 10), Pg = 1.10 x 10~ 4 ,and Oh = 3.14 x 10~ 4 . To begin, 
Fig. 10(a) exhibits an isolated, damped magnetic mode (labeled M), which, due to the singular nature of the 
limit Pm — > 0, is entirely absent from Fig. 10(b). Its magnetic energy E^jE — 0.28628 is the largest of the 
modes with Im(c) > — 1. Mode M, which is also present in channel problems at comparable {Re, H z , Pm, a}, 
has sufficiently large decay rate so as not to affect the validity of the inductionless approximation in critical 
Reynolds number calculations (see [19] and §5.4 ahead). The second contrastive feature is the presence of 
an unstable P mode ((U) < Re(c) < 1) in the Pm = 10~ 4 spectrum, when all of the modes of the induction- 
less problem are stable. At the same time, the decay rate —T = 0.020234a of the F mode (Im(c) > 1) in 
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Fig. 8. Eigenvalue spectra of film MHD problems at Re = 10 4 , Pm = 1.2, a = 1, and p u = Pb = 500. The external magnetic 
field is flow-normal (H x = 0), with H z = 14 (a) and H z = 100 (b). o and + markers respectively represent hydrodynamic and 
magnetic modes. Boldface markers correspond to unstable modes. 



Fig. 10(a) is about a factor of ten greater than its — r = —0.0035581a value in the inductionless limit. The 
magnetic energy of the Pi and F modes at Pm, = 10 -4 , respectively amounting to 0.21534S and 0.20240.E, 
is somewhat smaller than that of the M mode, but still more than an order of magnitude greater than the 
magnetic energy E b /E < 0.0064 of the remaining modes with Im(c) > —1. For these latter modes, the 
relative error in c of the inductionless approximation is less than 0.0058. 

The unstable P mode can be continuously traced to the F mode as Pm — ► 0, and likewise the F mode at 
Pm = 10~ 4 originates from the Pi mode in the inductionless limit. The relative change in c accumulated 
in the process is 0.027 and 0.021, respectively for Pi and F (in the sense of the Pm — 10~ 4 problem). 
This type of exchange of the modes' physical character, oftentimes accompanied by instabilities, is common 
in multiply diffusive systems [56] . Here, the velocity and magnetic-field perturbations play the role of two 
diffusive substances, respectively with diffusion constants i?e _1 and Rm~ . As discussed in more detail in [8], 
the free-surface is essential to this low- Pm instability, which does not occur in channel problems. 

The markedly different types of behavior we have so far encountered are a testament that in its full 
generality the free-surface MHD stability problem is a complex one. One of its major aspects that we have 
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Fig. 9. Same spectrum as in Fig. 8 but with H x = _ff z /tan(l°) « 5,792 
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Fig. 10. Spectra of film problems with flow-normal background magnetic field (H x = 0) at Re 
Oh = 3.14 X 10~ 4 , a = 0.01, Pm = 10~ 4 , p u = Pb = 500 (a), and the corresponding inductionlcss 
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not touched upon, and which we defer to future work, is the role of the induced magnetic field B on the 
instabilities, as well as the formation of multiple eigenvalue branches, at Pm = 0(1). While we do not 
present these calculations here, setting B to zero while keeping all other parameters fixed yields spectra that 
neither contain unstable magnetic modes, nor exhibit the multiple-branch structures. Using the approach 
employed in [8] for low- Pm fluids, it would be interesting to investigate the energy-transfer mechanisms 
associated with B, and the manner in which they contribute to the instability. 
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5.2. Convergence and Stability 

The issue of convergence and stability of spectral schemes is a very broad one, and can be approached from 
various angles. As a minimum, certain analytical criteria must be satisfied (e.g. §2.2 in [36]). That is, given a 
well-posed variational formulation of the problem at hand, the discrete solution must converge, under some 
suitable norm, to the exact one as the dimension N tends to infinity. Furthermore, the discretization error 
must be bounded by an TV-independent constant (stability). Among the relevant literature for eigenvalue 
problems (see [38] and references therein) of particular importance to us is the work of Melenk et al. [7], 
who showed that the Galerkin method used for what we call here non-MHD channel problems is spectrally 
convergent. Generalizing the results in [7] to MHD is of course an essential prerequisite if our schemes are to 
be deemed well-posed. In what follows, however, instead of pursuing that program we adopt a less rigorous 
approach and limit ourselves to more practical aspects of convergence and stability. That is, implicitly 
assuming that our schemes convergence in the analytical sense, we perform test calculations that aim to 
probe their behavior in actual computing environments, emphasizing on issues related to finite arithmetic 
precision. 

5.2.1. p- Convergence 

In shear-flow stability problems at large Reynolds numbers both truncation and roundoff errors come 
into play, and in certain cases addressing them leads to self-conflicting situations. On one hand, in order 
to resolve the small length scales that develop (the boundary and internal friction layers) it is necessary 
to work at large spectral orders (p > 500). Otherwise, the truncation error is significant. However, unless 
the basis polynomials are carefully chosen, the matrix representations of high-order differential operators 
(such as the D operator in the OS equation) become ill conditioned as p increases, causing a growth in 
roundoff error to the point where it exceeds truncation error. It is precisely here that a major strength of the 
scheme developed by KMS for plane Poiseuille flow, and extended here to free-surface MHD problems, lies: 
By working in polynomial bases constructed so as to reflect the order of the Sobolev spaces of the underlying 
continuous problem (see Remark 5), roundoff sensitivity essentially becomes independent of p. 

As a concrete illustration, we have experimented with an alternative implementation of our Galerkin 

[21 

schemes for non-MHD channel flow, where instead of the A„ polynomials prescribed in Table 1, the basis 
polynomials of V™" are Lagrange interpolants on LGL quadrature knots of order p + 1, suitably modified to 
meet the essential boundary conditions (2.3). Basis polynomials of this type, hereafter denoted by h n , are 
widely used in pseudospectral and spectral-element methods [36,57]. However, they lack the orthogonality 
properties appropriate to Hq. 

Remark 15. A prominent manifestation of non-orthogonality in the {h n } basis is matrix coefficient 
growth with p. We observed that the oo-norm of matrices with elements (D d2 (h n ), D^ 1 (/i m ))o,fi scales as 
pdi+d 2 j n con t r ast, all of the corresponding matrices T^ kdld ^ evaluated in the {An'}, {/i n } and bases 
(see Appendix A) have p-independent oo-norms. Recalling that the oo-norm of a matrix A is equal to 
max m ^2 n |^4 mn |, the latter is a direct consequence of the orthogonality properties of the Legendrc polyno- 
mials and the choice of normalization, which ensure that T^ kdld ^ (i) are banded, (ii) their bandwidths are 
p- independent, and (iii) apart from those corresponding to the nodal shape functions, the absolute values 
of the matrix coefficients cither remain constant or decrease going down the nonzero diagonals. 

The ill behaved stiffness and mass matrices arising in the Lagrange-interpolant basis lead to a rapid in- 
crease of the scheme's roundoff sensitivity with p. The resulting degradation in the accuracy of the computed 
eigenvalues is immediately obvious in Fig. 11, which shows the relative convergence of the least-stable eigen- 
value at {Re, a) — (10 4 , 1) as a function of the polynomial degree p, obtained via the {h n } and {Al?} bases. 
In both cases, convergence has been computed relative to a reference value obtained by means of the {A1T J } 
basis at high polynomial degree (p — 5,000). At small to moderate values of p the results are essentially 
identical, and clearly display the exponential decrease in truncation error typical of spectral methods. How- 
ever, in the case of the eigenvalue computed using the {h n } basis the exponential convergence trend halts 
abruptly at around p = 50, at which point the roundoff error caused by the ill conditioned stiffness and mass 
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Fig. 11. Spectral convergence of the least stable eigenvalue (the Ai eigenvalue in Table B.l) of non-MHD channel flow with 

the Poiseuille velocity profile at Re = 10 4 and a = 1. The solid and dotted lines respectively correspond to eigenvalues 

[21 

computed using the {A„ } basis (4.10a) with polynomial degree p, and the {h n } basis, consisting of Lagrange interpolants at 
the Legendre-Gauss-Lobatto quadrature knots of order p + 1, modified to satisfy the Hq boundary conditions. Convergence is 

evaluated relative to a reference value computed using the {a[?} basis at p = 5,000. 

matrices becomes dominant. As p further increases, the eigenvalue is seen to progressively diverge from the 
reference value until, around p = 400, the algorithm for the computation of the differentiation matrices (see 
Appendix C in [57]) becomes unstable and breaks down. In contrast, the eigenvalue computed using the 
{A„ } basis converges exponentially until close to machine precision, and although a small systematic trend 
can be observed for p > 10 3 , the calculation remains stable even at very large p. 

5.2.2. Effects of the Hartmann Profile 

Problems with the Hartmann velocity and magnetic-field profiles (2.12) differ from their counterparts with 
quadratic (or, more generally, polynomial) steady-state profiles in that the stiffness matrix K contains con- 
tributions of the form d£ e ff ^L m (£)£„(£), where is a real parameter. These exponentially-weighted 
inner products are nonzero for all (m,n), and, as a result, K is full. One immediate implication concerns 
memory and eigenvalue-computation costs, respectively scaling as N 2 and N 3 for a problem of dimension 
N = dim(V r [ Jv l). MHD problems are especially affected, since the spectral decompositions now have to be 
performed for both of the velocity and magnetic-field eigenfunctions, leading to 4-fold and 8-fold increases in 
size and complexity relative to inductionless or non-MHD problems. The non-sparsity of K also necessitates 
a re-evaluation of whether or not our schemes are roundoff stable at large spectral orders. For, our argument 
in Remark 15 that H-RTHoo is p-independent relied on the number of nonzero elements in each of its rows 
being fixed, which no longer applies in problems with exponential profiles. Yet, as Fig. 12 illustrates, in 
practice H-K^loo is to a very good approximation p- independent irrespective of the value of the Hartmann 
number, suggesting that our schemes are well-conditioned for the Hartmann family of steady-state profiles 
as well. Of course, ||-K"||oo does experience a growth with H z , but that growth is due to physical parameters 
only. 

In §4.2.2, we introduced two alternative ways of evaluating the U and independent terms in the stiffness 
matrix, one of which employs suitable quadrature rules [34] to compute the exponentially-weighted inner 
products exactly, while the other is based on approximate LGL quadrature at the precision level specified 
in (4.44). The eigenvalue calculations in Table B.3 have already hinted at a close agreement between the 
two methods in inductionless flow, which we now examine in more detail, using film MHD flow with oblique 
magnetic field as a more challenging example. We consider a problem with the same parameters as in 
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Fig. 12. Infinity norm of the stiffness matrix K of film MHD problems with the Hartmann velocity and magnetic-field profiles 
at Re = 10 4 , Pm = 1.2, Pg = 1.10 X lO" 4 , Oh = 3.14 X lO" 4 , a = 1, and H x = H*/tan(l°). In the left-hand panel, ||.K"||oo 
is plotted as a function of p = p u = Pb a t H z = 0.1, 10, 500. The right-hand panel shows 1 1 fC| | oo as a function of H z at fixed 
Pu = Pb = 200. The dashed line represents the power law H-K^loo oc H 2 J A . 

Fig. 9, and track the dependence of the computed eigenvalue of Mode 1 and Mode 31 (as usual, ordered 
in descending order of Re(7)) as p = p u = Pb is varied from 30 to 1,500. We calculate the eigenvalues 
using both exact quadrature (Eqs. (4.40) and (4.41)), and approximate quadrature (Eqs. (4.46)-(4.47)) 
at the smallest precision level consistent with (4.44). Fig. 13 demonstrates that the eigenvalues converge 
exponentially towards their reference values, computed at p = 2,500 via the exact-quadrature method, in a 
nearly identical manner, until limited by finite arithmetic precision. Convergence for Mode 31 is about an 
order of magnitude less than Mode 1, but in both cases the computed eigenvalues remain stable at large p. 
It therefore appears that a version of Banerjee and Osborn's theorem [37] that 2p — 1 quadrature precision 
is sufficient for convergence in elliptical eigenvalues problems also applies in OS-type problems. We remark 
that due to the aforementioned issues regarding storage and computation cost, we were not able to extend 
the calculation to as high values of p as we did in the non-MHD problem with the Poiscuillc velocity profile 
(Fig. 11). 

5.2.3. Non-Normality Issues 

Despite yielding stiffness and mass matrices that are 'optimally' conditioned with p, our choice of bases 
does comparatively little in addressing the second major source of roundoff error in our stability problems, 
which is due to the non-normality of the OS and induction operators (2.1). As already discussed in §1.1, at 
large Reynolds numbers the OS operator is highly non-normal, and, in consequence, its spectrum contains 
nearly linearly dependent eigenfunctions (with respect to the L 2 or energy inner products). According to 
Rcddy ct al. [33], expanding arbitrary functions of unit norm in terms of the OS eigenfunctions would require 
coefficients scaling as exp(7i?e 1 ^ 2 ) (for a = 1). At around Re ~ Ax 10 4 , the coefficients would be as large 
as 10 16 , indicating that in 64-bit arithmetic (15 significant digits) expansions of arbitrary functions would 
be severely affected by roundoff error. Similarly, one would expect the reverse operation of decomposing the 
OS eigenfunctions in a basis of polynomials to be also characterized by a sharp rise in roundoff sensitivity 
with Re. 

Consider, for example, the spectra in Fig. 14, which have been computed at Re = 4 x 10 4 and Re = 10 5 
with our Matlab code, working in 64-bit arithmetic. Instead of a well-defined intersection point between the 
A, P, and S branches, the numerically computed spectra exhibit a diamond-shaped structure of eigenvalues, 
whose area on the complex plane grows with Re. This type of spectral instability, which is entirely caused 
by finite-precision arithmetic, has come to be the hallmark of roundoff sensitivity due to non-normality of 
the OS operator [11,12,31,33]. As expected from the analysis in [33], the sensitivity increases steeply with 
the Reynolds number: Comparing the spectrum at Re = 4 x 10 4 with the corresponding one at Re = 3 x 10 4 



33 



10 ' 1 1 — 1 — 1 — 1 — 

2 

10 



Fig. 13. p-convergence in film MHD problems with Hartmann velocity and magnetic-field profiles and oblique steady-state 
magnetic field. The flow parameters are as in Fig. 9 and Table B.7, and the eigenvalues shown are for Modes 1, and 31. The 
solid and dashed lines respectively represent eigenvalues obtained via the exact and LGL quadrature schemes. Convergence is 
computed relative to a reference value at p = p u = p;, = 2,500 obtained using the exact-quadrature method. 




Re(c) Re(c) 

Fig. 14. Eigenvalue spectra of non- MHD film problems at Pg = 1.1 X 10~ 4 , Oh = 3.14 X 10~ 4 , a = 1, p u = 500, and Reynolds 
numbers Re = 4 X 10 4 (a) and 10 5 (b) 

(Fig. 4) reveals that it only takes a factor of 0.3 increase in Re for a noticeable diamond-shaped pattern to 
form (though a small diamond is already present in Fig. 4). 

In channel problems, the diamond-shaped pattern also emerges at around Re = Ax 10 4 , and at marginally 
smaller i?e = 3x 10 4 if the calculation is performed in the Lagrange interpolant basis of §5.2. In the latter 
case, decreasing Re to 2 x 10 4 is sufficient for the diamond to become virtually unnoticeable by eye, despite 
the basis being 'ill-conditioned'. A similar Reynolds number for the onset of the spectral instability {Re = 
2.7 x 10 4 in Fig. 2 of Ref. [12]) is reported by Dongarra et al. for their Chebyshev tau scheme. Therefore, 
in these examples the accuracy of the computed eigenvalues and eigenvectors appears to be limited by the 
physical parameters of the problem, in accordance with the estimates of Reddy et al. [33] , rather than the 
details of the numerical scheme. If that is the case, then, as noted by Dongarra et al. [12], the only way of 
addressing the non-normality issue would be to increase numerical precision. Those authors have observed 
that working in 128- bit arithmetic does indeed remove the diamond-shaped pattern from the numerical 
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Fig. 15. Eigenvalue spectra of film MHD flow with the Poiseuille velocity profile and zero steady-state magnetic field 
(H x = H z = 0) at magnetic Prandtl numbers Pm = 5 (a) and 10 (b). The remaining parameters are equal to those in 
Fig. 7. As Pm increases, the magnetic modes (marked with + markers) develop the diamond-shaped pattern characteristic to 
roundoff errors caused by non-normality of the stability operators. The hydrodynamic modes (represented by o markers) are 
accurately computed, as they are decoupled from the magnetic ones (and do not depend on Pm). 

spectra. Unfortunately, we have not been able to verify this for our Galerkin schemes, as our code was 
written in Matlab, which does not natively support extended-precision floating-point numbers. However, 
there is no reason to believe that increasing the number of significant digits would not alleviate the spectral 
instability in our schemes as well. 

Turning now to MHD, at a given value of the Reynolds number the effects of non- normality may be more or 
less severe compared to the hydrodynamic case, depending on the remaining parameters of the problem (Pm, 
H x , H z ). The general rule of thumb is, however, that whenever the spectrum contains branch- intersection 
points, the highly non-orthogonal modes close to them will at some point experience the spectral instability 
if Re and/or Rm arc increased. The examples in Fig. 15 illustrate that in problems with zero background 
magnetic field the magnetic modes are the first to develop the diamond-shaped pattern if Pm is greater 
than unity. Moreover, Fig. 16 shows that if Re is increased in the film MHD problem in Fig. 9 four branch 
intersection points are formed, all of which are affected by roundoff errors at Re = O(10 5 ). On the other 
hand, in the Pm < 10~ 5 regime relevant to terrestrial fluids, the gradual disappearance of the three-branch 
structure with increasing H z (see §5.1.2) results to smaller regions on the complex plane being dominated 
by inaccurately computed eigenvalues. 



5.3. Consistency Calculations 

The energy-balance relation (2.17) forms the basis of the following consistency check for film MHD prob- 
lems: First, solve the matrix eigenproblem (4.3) to obtain 7 and the discrete representation v of (u,b,a). 
Using v, the corresponding definition of the basis functions tjj m (Def. 6), and Eqs. (2.16), compute the 
quantity r :— Jr + 7m + r,j + P v + r v + P av + r arj . Then, according to (2.17), the relative difference 
between r = Re(7) and 7, given by e := \(t — r)/T\, should be small, ideally close to machine precision. 
We note that the presence of the second derivatives of b in (2.16e) and (2.16g), which cannot be defined 
weakly for b E 77 1 (J7), necessitates that for the purposes of this calculation (u, b, a) is restricted to the 
strong solution space T>jc- Of course, in our polynomial subspaces of 77 1 square integrability of second (and 
higher) derivatives is in principle not an issue. However, as we discuss below, practical repercussions in 
evaluating expressions like (2.16e) and (2.16g) are nonetheless present, since in the {/i„} basis, which has 
been constructed so as to reflect 77 1 regularity, the matrix representations of scsquilincar forms involving 

second weak derivatives are not stable with p. Specifically, as can be checked either numerically or from the 

- 2 - 2 

properties of the Legendre polynomials, matrix coefficients of the form (D fj, n , D /U m ) q grow like p 2 , while 
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Fig. 16. Eigenvalue spectra of film MHD flow with oblique steady-state magnetic field at Reynolds numbers Re = 5 X 10 4 
(a) and Re = 10 s (b). Apart from Re, all other flow and spectral parameters are equal to those in Fig. 9. In these plots no 
distinction is made between hydrodynamic and magnetic modes. 

" 2 

boundary terms D /Xn(l) scale as p 5 ^ 2 . 

Fig. 17 shows the details of such a calculation for film MHD problems at Pm = 1.2, with external magnetic 
field oriented at 1° relative to the streamwise direction, and flow-normal Hartmann number H z E [0.1, 100] 
(correspondingly, the streamwise Hartmann number H x ranges from approximately 5.73 to 5,730). As H z is 
varied, a single mode is continuously tracked, which corresponds to Mode M 2 at H z — 100 (sec Table B.7 
and Fig. 9). That mode is stable for sufficiently weak magnetic fields, but as H z increases it undergoes 
an instability in which the dominant power input is Maxwell stress (cf. the instabilities in non-MHD and 
low- Pm flows caused by positive Reynolds stress). At the same time, the energy E (2.15) changes from being 
predominantly magnetic to a nearly equal mix of magnetic and free-surface energies (at H z <~ 10 the energy 
is also seen to have a significant kinetic contribution). For all values of the Hartmann number considered, 
the error e remains small (e < 10~ 6 ), but displays a trend with H z that mirrors r aj} . We attribute this 
behavior to roundoff error in r ar/ due to that term's dependence on D 2 6(0). In fact, the reason that we 
chose to examine Mode M 2 , rather than, say, Mode M 1; is that at sufficiently large Hartmann numbers 
the magnetic and surface energies of that mode are both appreciable, making it particularly susceptible to 
errors associated with r ar) . Indeed, as the dotted line in the lower- left panel in Fig. 17 shows, decreasing p 
from 200 to 100 results to a noticeable change in e, which diminishes roughly by an order of magnitude. On 
the other hand, modes with small \r ar) \, are comparatively unaffected by the choice of p (e.g. for Mode Mi 
e ~ 10~ 10 for both p = 100 and p — 200, and for all H z < 100). It therefore appears that e is dominated by 
roundoff error in r ari , rather than some inconsistency in our numerical scheme and/or its implementation. 

As a further consistency check we have compared growth rates for the two most unstable modes of the 
results of Table B.2 at Re = 3 x 10 4 (see also Fig. 4(b)) with the results of free-surface flows computed using a 
fully nonlinear Navier-Stokes solver. The Navier-Stokes code is based on the arbitrary Lagrangian-Eulerian 
spectral element code developed by Ho [15], Ronquist [58], and Fischer [59]. For a = 1, the (nominal) 
computational domain was taken as fl — [0, 2tt] x [—1,0], which was tessellated with a 6 x 10 array of 
spectral elements. A uniform element distribution was used in the streamwise direction while a stretched 
distribution was used in the wall-normal direction. Near the wall, an element thickness of Az — 0.005 was 
used to resolve the boundary layer of the unstable eigenmodes. The polynomial order within each element 
was N = 13 and third-order timestepping was used with At = 0.00125. The initial conditions corresponded 
to the base flow plus S := 10~ 5 times the velocity eigenmode associated with mode k, with k = 1 or 2, 
i.e. the most and second-most unstable eigenmode for these particular flow conditions. The domain was 
stretched in the z direction to accommodate the 0(<5) surface displacement using transfinitc interpolation 
[60]. The eigenmodes, which are defined only on z = [—1,0], were mapped onto the nominal domain then 
displaced along with the mesh. The base flow was defined as U(z) = 1 — z 2 over the deformed mesh. 
Mean growth rates were computed by monitoring the L 2 -norm of the wall-normal velocity and defining 
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Fig. 17. Energy balance for film MHD flow at Re = 10 4 , Pm = 1.2, Pg = 1.1 X 1CT 4 , Oft = 3.14 X 1CT 4 , a = 1, and 
Pu = Pb = 200. The flow-normal Hartmann number H z ranges from 0.1 to 100, with the streamwisc Hartmann number given 
by H x = H z / tan(l°). The curves track the behavior of a single mode as a function of H z , which at H z = 100 is Mode M2 
(see Fig. 9 and Table B.7). The graphs in the right-hand panels show the energies (2.15), normalized so that E = 1, and the 
power-transfer terms (2.16). The solid (dotted) portions of the curves in the logarithmic plots correspond to positive (negative) 
values. The left-hand panels display the growth rate r and phase velocity C, as well as the error e. The latter has also been 
evaluated at p u = Pb = 100, and plotted as a dotted line, in order to illustrate the roundoff sensitivity in r a ri- 



r(t) := ln(||w(i)|| 2 /||w(0)||2). The error is again denned as (-T(i) — r)/T, where r is computed using the 
linearized code. Aside from some initial transients, the error over t = [5,100] was less than 10~ 5 for Mode 1 
(r = 0.007984943826437) and less than 5~ 3 for Mode 2 (T = 0.000052447145102). These results provide 
independent confirmation of both the linear-stability and spectral-element based ALE codes. 
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Table 2 

Hartmann- number dependence of the critical Reynolds number Rc, wavenumbcr ct c , and phase velocity C c for the even mode 
in channel problems, and the hard and soft modes in film problems. The critical parameters have been computed in the 
inductionless limit and Pm = 10~ 4 , for H x = and H z £ [0, 100]. N u = JVj, is the dimension of the velocity and magnetic-field 
solution spaces used in the calculations. The underlined digits in the results for the even channel mode indicate discrepancy 
from the corresponding calculations in Tables 1 and 3 of Ref. [19]. The soft mode's critical parameters in the inductionless limit 
were evaluated via Eqs. (5.1). 

Inductionless Pm = 10~ 4 
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5.4. Critical-Parameter Calculations 

Our final set of calculations pertains to the critical parameters for the onset of instability in channel and 
film problems with flow-normal background magnetic field and Hartmann steady-state profiles. In channel 
problems, we seek the minimum (critical) Reynolds number Rc, and the corresponding wavenumber a c , 
for which the spectrum contains unstable modes, keeping the Hartmann number H z and, where applicable, 
the magnetic Prandtl number Pm fixed. In film problems, we also constrain the Ohnesorge number and 
the gravitational Prandtl number, setting Pg = 1.10 x 10 -4 and Oh = 3.14 x 10~ 4 . As stated in §2.2, 
neither of H z , Pm, Pg and Oh depend on the velocity scale of the flow. Therefore, our formulation of the 
critical-parameter calculation leads, among other possibilities, to a determination of the minimum steady- 
state velocity at which the flow becomes linearly unstable when all of its remaining properties are held 
fixed. Other choices of parameter constraints may be appropriate depending on the particular application. 
In all cases, however, Rc, a c , and the corresponding modal phase velocity C c , can be obtained by solving 
a minimization problem for Re, constrained by the cigcnproblcm (4.3) and the normalization ||u||2 = 1- 
The numerical results presented in Tables 2 and 3 were obtained in that manner, using Matlab's fmincon 
optimization solver to carry out the computations. 

We begin from channel and film problems in the inductionless limit, critical parameters of which are listed 
for H z G [0, 100] in the left-hand portion of Table 2. In the channel case, the critical mode is always of even 
symmetry and lies in the A branch of the spectrum (C c < (U)). As indicated by the underlined digits in the 
computed parameter values, our calculations are in excellent agreement with those by Takashima [19]. Even 
though channel problems also exhibit an odd unstable mode, its critical Reynolds number always that of 
the even one [54], and therefore we do not consider it here. On the other hand, in film problems cither the 
soft or the hard mode, respectively characterized by C < (U) and C > 1 (see §5.1), can have the smallest 
critical Reynolds number, depending on the Hartmann number and the free-surface parameters. The critical 



38 



Table 3 

Critical Reynolds number Rc, wavenumber ct c , and phase velocity C c for channel and film MHD problems at (H x , H z ) = (0, 10) 
and Pm € [10~ 8 , 10~ 4 ]. All calculations were performed at dimension N u = Nf, = 300. The underlined digits in the channel- 
problem calculations differ from the corresponding ones in Table 3 of Ref. [19] . 

Channel Film 



Pm 




Rc 


Olc 


C c 




Rc 




a c 




Cc 




1.0E 


-08 


4.3981789E + 05 


1.739135E + 00 


1.547887E 


- 01 


4.3978679E - 


-05 


1.739235E 


+ 00 


1.547884E 


- 01 


1.0E 


-07 


4.3981547E + 05 


1.739135E + 00 


1.547889E 


- 01 


4.3978440E + 05 


1.739234E 


+ 00 


1.547886E 


- 01 


1.0E 


-06 


4.3979162E + 05 


1.739128E + 00 


1.547912E 


-01 


4.3976064E - 


-05 


1.739228E 


+ 00 


1.547909E 


- 01 


1.0E 


- 05 


4.3958738E + 05 


1.739141E + 00 


1.548125E 


- 01 


2.1862147E- 


-05 


1.947639E 


- 03 


1.004972E 


+ 00 


1.0E 


- 04 


4.3861946E + 05 


1.739024E + 00 


1.549340E 


- 01 


1.0019578E- 


-05 


3.579587E 


- 03 


1.015528E 


+ 00 


1.0E 


- 03 


4.2969213E + 05 


1.739870E + 00 


1.559585E 


- 01 


6.1108706E- 


-04 


2.914764E 


- 03 


1.020771E 


+ 00 


1.0E 


- 02 


4.8282141E + 04 


4.894029E - 03 


8.973103E 


- 01 


4.6367226E - 


-04 


9.120983E 


- 04 


1.089575E + 00 


1.0E 


- 01 


6.8382770E + 02 


2.788195E - 01 


8.899146E 


- 01 


1.1205597E- 


-03 


1.869298E 


- 01 


8.803851E 


- 01 



wavenumber of the soft mode is zero, and in light of this it is possible to derive the closed-form expressions 

Rc = 23/ 2 sinh( gz /2)( g -tanh( gz ))V 2 C c = 1 + sech(^), (5.1) 

Pg(H z coth(i? z /2) sech 3 (iJ z )(2ff z (2 + cosh(2i? z )) - 3 sinh(2 J ff z ))) 1 /2 

using regular perturbation theory about a = [8]. In the H z — > limit, Eqs. (5.1) reduce to Rc = (5/8) 1 / 2 / Pg 
and C c = 2. At H z > the critical Reynolds number grows exponentially, implying that for all but small 
Hartmann numbers the onset of instability in inductionless film problems is governed by the hard mode. 
Direct numerical calculations of the soft mode's critical parameters rapidly become intractable, but wc 
checked that the linear-stability code yields results of O(10~ 8 ) accuracy for H z < 10. As for the hard mode, 
the inductionless results in Table 2 show that apart from small Hartmann numbers, where gravity and 
surface tension are more important than the magnetic field, its critical parameters are very close to those 
of the channel mode, suggesting that for sufficiently strong magnetic fields the free surface only plays a 
minor role in the hard instability. The critical wavenumber of the channel and hard modes increases with 
H z (i.e. shorter wavelengths become unstable first), which is consistent with the fact that the decreasing 
Hartmann-layer thickness is the principal contributing factor in the instability suppression [54] . 

We now examine how nonzero magnetic Prandtl numbers modify the picture in the inductionless limit. 
As can be checked from Table 2, the error in Rc incurred by making the inductionless approximation is 
less than 4 x 10 -3 for the even channel mode over all Hartmann numbers probed. Moreover, according to 
Table 3, at H z — 10 the critical Reynolds number of the even channel mode decreases by a factor smaller 
than 0.003 when Pm is increased from 10 -8 to 10~ 4 . These calculations are in very good agreement with 
the corresponding ones by Takashima, and are illustrative of the weak dependence of critical parameters 
of the channel mode on Pm < 10~ 4 for all H z < 200 [19]. As Pm grows above O(10~ 4 ), the accuracy of 
the inductionless approximation progressively deteriorates, until the critical mode undergoes a bifurcation 
to a magnetic mode (i.e. a singular mode in the limit Pm — > 0) of odd symmetry, manifested by the sharp 
decrease in a c at Pm = 10~ 2 . 

Turning to film problems, Table 2 demonstrates that, as with the even channel mode, at small Pm the 
inductionless approximation yields accurate results for the critical parameters of the hard mode. On the 
other hand, the data clearly show that a small, but nonzero, Pm affects profoundly the critical parameters 
of the soft mode. In particular, the previously encountered exponential growth of Rc with H z becomes 
suppressed to the point that it now trails the hard mode's critical Reynolds number by a wide margin. In 
the right-hand portion of Table 3 the hard mode (C c < 1) is seen to govern the onset of instability for 
Pm < 10~ 6 , with the soft one, characterized by C c > 1, taking over at larger magnetic Prandtl numbers. 
Even though no further bifurcations occur for Pm e [10~ 6 ,10 -2 ], the soft mode in itself appears to be 
sensitive to Pm. In the context of large-wavelength perturbation theory, this behavior can be traced to 
the coefficient 71 in the asymptotic expansion 7 = cryi + a 2 72 + 0(a 3 ), which vanishes as Pm — > 0. In 
consequence, the equation 71 + aj2 = 0, used to determine Rc, becomes singularly perturbed, resulting 
to the observed sensitivity. From a physics standpoint, as already mentioned in §5.1.3 and discussed in 
further detail in [8], the problems with nonzero Pm are susceptible to doubly diffusive effects, giving rise 
to instabilities not present in the inductionless limit. In total, over the interval 10~ 8 < Pm < 10 -4 , which 



39 



roughly coincides with the Pm values of terrestrial incompressible fluids, the critical Reynolds number of 
the examined film problems decreases by almost a factor of five. 

6. Conclusions 

In this paper we have presented a spectral Galerkin method for linear-stability problems in free-surface 
MHD. The method is essentially an extension of the scheme developed by Kirchncr [6], and Melenk, Kirchncr 
and Schwab [7] to solve the Orr-Sommcrfeld (OS) equation for plane Poiscuillc flow. Besides free-surface 
MHD problems, which we refer to as film MHD problems (Dcf. 1), our scheme provides a unified framework 
to solve MHD stability problems with fixed boundaries — the so-called channel MHD problems (Def. 2) — and 
their simplified versions at vanishing magnetic Prandtl number Pm, which we refer to as inductionless film 
and channel problems (Defs. 3 and 4). We studied problems with either the Poiseuille velocity profile, or 
the Hartmann velocity and magnetic-field profiles, both of which are physically motivated. However, our 
schemes are applicable to arbitrary analytic steady-state profiles. In all cases, the Galerkin discretization 
results to the matrix generalized eigenvalue problem Kv = jMv, where K and M are respectively the 
stiffness and mass matrices, 7 is the complex growth rate, and v is a column vector containing the problem's 
degrees of freedom. We detected no spurious eigenvalues, a fact which we attribute to the non-singularity of 
M in all bases. 

We discretized the solution spaces for the velocity and magnetic-field cigenfunctions using Lcgcndre in- 
ternal shape functions and nodal shape functions, chosen according to the Sobolev spaces of the continuous 
problems. Separating the basis polynomials into internal and nodal ones facilitates the natural (weak) im- 
positition of the boundary conditions for free-surface MHD, namely the stress and kinematic conditions at 
the free surface, and the Robin-type insulating boundary conditions for the magnetic field. The orthogonality 
properties of the bases guarantee that roundoff error is independent of the spectral order p, allowing one 
to work at the large spectral orders (p > 500) required to resolve the small length scales present at high 
Reynolds numbers Re. Moreover, in problems with polynomial velocity and magnetic- field profiles, K and 
M are sparse, and iterative solvers can be used to compute 7 and v efficiently. 

The optimal conditioning of our schemes with respect to p alleviates only marginally their roundoff 
sensitivity due to non-normality of the stability operators. At around Re — 4x 10 4 we observed the formation 
of the characteristic diamond shaped pattern on the complex eigenvalue plane caused by lack of sufficient 
precision in 64-bit arithmetic. An alternative discretization, performed in terms of Lagrange interpolation 
polynomials, was found to give rise to the pattern at only slightly smaller Reynolds numbers (Re = 3x 10 4 , 
which is close to the value reported in [12] for a Chebyshev tau scheme), despite the ill conditioning of the 
Lagrange interpolant basis. Roundoff errors associated to non-normality therefore appear to be governed 
by physical parameters, rather than the details of the discretization scheme. Working in extended-precision 
(e.g. 128-bit) arithmetic is probably the only way to address this type of error, but, at the time of writing, 
that option could not be implemented with our Matlab code. 

We described two ways of addressing the presence of exponentially weighted sesquilinear forms in prob- 
lems with Hartmann steady-state profiles. In the first approach, the forms are evaluated without incurring 
quadrature error by means of the algorithm developed by Mach [34] to compute Gauss quadrature knots 
and weights for exponential weight functions on a finite interval. The second approach involves replacing the 
forms by approximate ones derived from Legendre-Gauss-Lobatto (LGL) quadrature rules at the 2p—l pre- 
cision level. The latter has been established by Banerjee and Osborn [37] as sufficient to guarantee stability 
and convergence in elliptical eigenvalue problems, but, to our knowledge, no corresponding bound exists for 
OS problems. We found that eigenvalues computed via the LGL method agree to within roundoff error with 
the corresponding ones obtained using exact quadrature, indicating that a version of Banerjee and Osborn's 
theorem should also be applicable in eigenvalue problems of the OS type. 

As an independent consistency check, we compared modal growth rates in non-MHD free-surface flow to 
energy growth rates in fully nonlinear simulations. At Re = 3 x 10 4 and wavenumbcr a=lwe found that 
the error over 100 convective times is less than 10~ 5 and 5~ 3 , respectively for the first and second least 
stable modes. We also compared modal growth rates in problems with oblique external magnetic field to the 
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corresponding ones derived from an energy conservation law for free-surface MHD. Here the error was found 
to be less than 1CP 6 , with its largest portion attributed to roundoff sensitivity in the calculation of one of 
the energy terms, rather than inconsistencies in the numerical scheme. In channel problems, we found that 
our results for the critical Reynolds number, wavenumber, and phase velocity for Hartmann flow agree very 
well with the corresponding ones by Takashima [19], obtained using a Chebyshev tau method. 

At the magnetic Prandtl number regime of terrestrial fluids (Pm < 10 ) and for H z > 5, the critical 
parameters of the hard instability mode in film flow were found to be close to those of the even critical 
mode in channel flow. Increasing Pm from 10 -8 to 10~ 4 at H z = 10 resulted to a mild, O(10~ 3 ), decrease 
of the critical Reynolds number Rc for the hard and channel modes, but Rc dropped by more than a factor 
of four for the soft (surface) instability mode in film flow. A surface- wave instability at small Pm, but 
absent in the inductionless limit, was also observed in the spectra of film MHD problems at aRe = 10 4 , 
H z — 10 and Pm = 10~ 4 . These results are indicative of the important role played by the working fluid's 
magnetic Prandtl number in the stability of industrial and laboratory free-surface flows. In test problems at 
Pm = 1.2 we observed that increasing H z from zero to 100 leads to the formation of multiple branches on the 
complex-eigenvalue plane. Unlike problems at small magnetic Prandtl numbers, the spectra at Pm = 0(1) 
contain unstable magnetic modes, two of which were recorded in a film MHD problem with oblique external 
magnetic field. 

Before closing, we note a number of directions for future work. On the analytical side, it would be highly 
desirable to extend the convergence analysis of Melenk et al. [7] to free-surface MHD. Even though the 
calculations presented in §5 provide strong numerical evidence that our schemes are stable and convergent, 
their wcll-poscdness cannot be settled without a rigorous analytical backing. Similarly, our proposed method 
in §4.2.2 of approximating weighted sesquilincar forms using LGL quadrature requires an adaptation of 
Banerjee and Osborn's [37] work to OS-type eigenvalue problems. A further analytical objective would be 
to generalize the criterion of scale resolution [7] , which provides an estimate of the minimum spectral order 
required to achieve convergence at a given Re in non-MHD channel flow. Of course, any such criterion 
would have to be numerically tested. On the physics side, our discussion in §5.1 and §5.4, which is mostly 
phenomenological, should be supplemented by a study of the operating physical mechanisms. In [8], we 
pursue such a study at the low-Pm regime, but that should be extended to cover Pm — 0(1) flows, which 
have been conjectured [5] to be relevant in certain astrophysical accretion phenomena. 
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Appendix A. Matrix Representations of the Schemes' Forms and Maps 

In this Appendix we provide expressions for the matrix representations of the sesquilinear forms and maps 
used in the main text. In §A.1-§A.3 we consider the T matrices, defined in (4.22) and (4.31), whose elements 
can be stably evaluated in closed form by means of the orthogonality properties of the Legendre polynomials 
(4.7). We then describe, in §A.4, how Mach's quadrature scheme [34] can be used to evaluate the matrices 
S (4.39) and C. 

A.l. The Matrices T l * d r ld2] 

We evaluate the matrices T 1 ^ 1 ^ e M. NxN (4.22a) listed in Table A.l. In the rightmost column of that 
table stands for the main diagonal and m (—to) represents the m-th upper (lower) diagonal. Also, in 
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Table A.l 



Properties 


of the matrices T'l^ 1 ^ 2 ^ 








r 




Symmetry 


Bandwidth 


Nonzero Diagonals 








S 








1 





c 
D 


2 


0, ±2 


1 


10 


s 


3 


±1, ±3 


1 


2 


s 


4 


0, ±2, ±4 


2 





s 


4 


0, ±2, ±4 


2 


10 


A 


3 


±1, ±3 


2 


10 


s 


5 


±1, ±3, ±5 


2 


110 


N/A 


4 


0, ±2, ±4 


2 


2 


S 


6 


0, ±2, ±4, ±6 



the third column from the left, S and A respectively identify symmetric and antisymmetric matrices. For 
each k and r, one only needs to evaluate the cases (di,d,2) — (0,0) and (1,0), since the results for the 
remaining values of d\ and di < r follow by making use of the hierarchical relation (4.23) and the property 



-\kd1d2] 



( rp[kd,2dl 



. Some of the results below can also be found in the paper by Kirchncr [6]. However, 

since that reference contains a number of typographical errors, and for the sake of completeness, we have 
opted to reproduce them here. 

Working down the rows of Table A.l, our first result, which has already been stated in (4.11), is simply 
r'°o°' = In- Next, we consider the nonzero elements in the main and upper diagonals of the matrices with 
r = 1, all of which are symmetric. These are 

1 



-,[000] 
H}, 



m = n — 2 : — 



v /(2n-3)(2n - l)V2n + l 
1/1 1 



m 



n : 



+ 



2n + 1 V 2n - 1 2n + 3 



(A.l) 



^[100] 

1 *l 



ran 

m = n — 3 : 
m = n — 1 :- 



V2n - 5(2n - 3)(2n - l)v / 2n+T 
1 / n- 1 



n 



+ 



n+l 



and 



^[200] 

if 1 

J ran 

m = n — 4 : — 
m = n — 2 : 



y/(2n - l)(2n + 1) V ( 2n - l )( 2n - 3) (2n-l)(2n+l) (2n+l)(2n + 3) 



(n - l)(n-2) 



(A.2) 



V2n- 7(2n - 5)(2n - 3)(2n - l)V2n + l 
1 / (n-2)(n-l) 



2(n- l) 2 

v /(2n-3)(2n+l) V(2n - 5)(2n - 3)(2n - 1) (2n-l)(2n+l) V ( 2 ^ - 3) 
n(n + 1) 



+ 1 



+ 



m = n: 



(2n- l)(2n+l)(2n + 3) y 
1 / 1 /2(n-l) 2 



(2n+l) V(2n-l)(2n + l) V (2n - 3) 
1 /2(n+l) 2 



2n(n + 1) 



(2n- l)(2n+ l)(2n + 3) 



+ 



(2n + 3)(2n + 5) V 2n + 1 



+ 1 



(A.3) 
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Among the r = 2 matrices, T^°\ T^ ' and T^ ' are symmetric, t'°2°' is antisymmetric, and T^ ' has 

no symmetry property. In Ref. [6], the matrix corresponding to T^ ' is denoted by T* . An expression for Ti 



is provided in that paper's Appendix B.6, but contains typographical errors. In Eq. (A. 7) ahead we indicate 
the erroneous terms, and also an error in the second upper diagonal (which we have been unable to trace to 
individual terms), by underlines. The nonzero elements of the symmetric and antisymmetric matrices, again 
in their and main and upper diagonals, are 



[000] 

-2 




mn 

m = n — 4 : 



V2n - 5(2n - 3)(2n - l)(2n + l) v / 2lT+3' 

-1/1 
m = n-2 : - - — — — — + — — — + 



V2n - l(2n + l)y/2n + S \{2n + 3)(2n + 5) (2n+l)(2n + 3) (2n-l)(2n+l) 

+ TTToZ os ) > ( A -4) 



(2n- l)(2n- 3) 

1/1 1 
m = n :- — — — — ■ , Nn/ „ — — — ■ + — — — — — + 



2n + 3\ (2n + 5) 2 (2n + 7) (2n + 3)(2n + 5) 2 (2n + l)(2n + 3)(2n + 5) 
1 1 

+ 



(2n+ l) 2 (2n + 3) (2n + l) 2 (2n - 1) 



^[010] 
- J mn 

m = n — 3 : 
m = n — 1 : 



1 



V2n- 3(2n - l)(2n + l)V2n + 3' 
1 / 1 



+ 



+ 



y/(2n + l){2n + 3) V ( 2 ™ - l)(2n + 1) (2n+l)(2n + 3) (2n + 3)(2n + 5) 



(A.5) 



-,[100] 

h„ 2 



m = n — 5 
m = n — 3 : 



n — 1 



V2n - 7(2n - 5)(2n - 3)(2n - l)(2n + l)V2n + 3' 
-1 / n-1 



+ 



n-1 



V2n - 3(2n - l)(2n + l)v^n + 3 V(2n - 5)(2n - 3) (2n - 3)(2n - 1) 
n n + 1 n + 1 



m = n — 1 : 



(2n-l)(2n+l) (2n+l)(2n + 3) (2n + 3)(2n + 5) / ' 
1 / n-1 2n 



(A.6) 



v /(2n+l)(2n + 3) V (2n - 3)(2n - l) 2 (2n + 1) (2n - l) 2 (2n + l)(2n + 3) 
4(n+l) 2(n + 2) 



(2n - l)(2n + l)(2n + 3)(2n + 5) (2n + l)(2n + 3)(2n + 5) 2 
n + 3 



(2n + 3)(2n + 5) 2 (2n + 7) / ' 



-,[200] 
Hi 
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m = n — 6 : 



m = n — 4 : 



(n-2)(n-l) 

V2n - 9(2n - 7)(2n - 5)(2ra - 3)(2n - l)(2ra + l)y/2n + 3' 

n(n + 1) 



V2n - 5(2n - 3)(2n - l)(2n + l)V2n + 3 V (2n + 3)(2n + 5) (2n+l)(2n + 3) 



n(n + 1) 



+ 



(n-1) 2 



(n-2)(n-l) (n-2)(n-l) 



m = n — 2 :- 



(2n-l)(2n+l) (2n-3)(2n-l) (2n-5)(2n-3) (2n - 7)(2n - 5) / ' 

(n + 2)(n + 3) 2 



2(n+l) 2 

V^n - l(2n + 1) V2n + 3 V ( 2 ™ + 3)(2n + 5) 2 (2n + 7) (2n + l)(2n + 5) 2 V 2n + 1 



+ 1 



+ - 



n(n + 1) 



+ 



+ 



(2n- l)(2n + 3) V(2n+l)(2n + 5) (2n-3)(2n + 5) (2n+l) 2 



+ 



1 



1 



m = n : 



2(n-l) 2 

(2n + l)(2n -3)7 (2n - 3)(2n + l) 2 V 2 ™ - 3 
(n-l)(n-2) 



(n + 3) 2 



+ 1 



(2n-5)(2n-3) 2 (2n-l)y ' 
(n + 4) 2 



2(n + 3)(n + 2) 



2n + 3 V (2n + 5) 2 (2n + 7) 2 (2n + 9) (2n + 5) 3 (2n + 7) 2 (2n + 3)(2n + 5) 3 (2n + 7) 
2(n + 3)(n + 2) (n + 2) 2 (n + l) 2 



+ 



+ 



(2n+ l)(2n + 3)(2n + 5) 2 (2n + 7) (2n + 3) 2 (2n + 5) 3 (2n + l)(2n + 3) 2 (2n + 5) 2 

2(n + 2) 2 2(n + l) 2 

(2n + l)(2n + 3) 2 (2n + 5) 2 + (2n + l) 2 (2n + 3) 2 (2n + 5) 
2n(n + l) 



+ 



(n + 2Y 



+ 



(n + 1) 2 



(2n- l)(2n+l) 2 (2n + 3)(2n + 5) (2n + l) 2 (2n + 3) 2 (2n + 5) (2n + l) 3 (2n + 3) 2 



2n(n + 1) 



(n — l) 2 



(2n- l)(2n+ l) 3 (2n + 3) (2n - l) 2 (2n + l) 3 (2n - 3)(2n - l) 2 (2n + l) 2 



Moreover, the nonzero elements of T^2°' are given by 



(A.7) 



-,[110] 
H 2 



m = n — 4 : 
m = n — 2 : 



n- 1 



v/(2n - 5)(2n - 3)(2n - l)(2n + l)y/(2n + 3) ' 

n — 1 



V(2n- l)(2n + 3) V (2n - 3)(2n - l)(2n + 1) (2n - l)(2n + l) 2 
(n + 1) (n + 1) 



(2n+ l) 2 (2n + 3) (2n + l)(2n + 3)(2n + 5) 



m = n :- 



n+1 



n + 1 



2n + 3\ (2n- l)(2n+ l) 2 (2n + l) 2 (2n + 3) (2n + l)(2n + 3)(2n + 5) 



n + 2 



n + 2 



+ 



n + 3 



(2n+l)(2n + 3)(2n + 5) (2n + 3)(2n + 5) 2 (2n + 5) 2 (2n + 7) 



(A.8) 



m = n + 2 : 



1 



n + 2 



+ 



n + 2 



m = n + 4 : 



^(2n + 3)(2n + 7) V ( 2 ™ + !)( 2 ™ + 3)( 2 n + 5) (2n + 3)(2n + 5) 2 
n+3 n+4 
~(2n+5) 2 (2n + 7) + (2n + 5)(2n + 7)(2n + 9) 
(n + 4) 



V2n + 3(2n + 5)(2n + 7)(2n + 9)V 2n + 11 ' 
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A.2. The Matrices T [ * dld2] and T [ ^ ld2] 

We now compute the matrices listed in Table A. 2. In light of Remark 9, we explicitly consider only the 
elements in their first two rows and columns with indices no greater than the spectral leakage /. The remaining 
elements can be deduced from the results in §A.l. All of the required T^f ld2 ^ matrices are symmetric. The 
nonzero elements in their first two rows are given by 



^,[000] 
1 H 1 



/ 2 1 1 1 

3 3 V6 3\/l0 
12 1 1 



\3 3 V6 3VW. 



^,[200] 

1 m 



^,[100] 
1 H 1 



-- — 1 

3 5V6 3-y/lO 15 



1 



( 4 


1 


1 1 


15 


15 


5\/6 7-/L0 


1 


4 


1 1 


V 15 


15 


5V6 7 VlO 



15 105 

V2 
105 J 



15 



1 1 

V 3 _ 5V6 _ 3VT0 "15"/ 



1 m 




(A.9a) 



(A.9b) 



where m < 2 and, in each case, n <l. The only-non symmetric y^i^] ma ^- r j ces are y^ ! j anc j T^ '- The 
nonzero elements in their first two rows are 



-,[010] 
H? 



/ 1 1 3 
2 ~5 7^10 







1 



105 V2 
1 







, 2 

1 Q A/5 1 

\5 ~2T 15V14 105V2 / 



(A.10) 



and 



-,[110] 



70 35 



1 



1 



\ 21 105 21\/T0 
where n < I. Moreover, for m < I and n < 2 we have 

/ 1 



15V14 




105V22 
1 



V2 

315 105V22 J 



-,[110] 

H? 



45 



-,[010] 







-,[010] 



and 



nm 
T 



315 



3V10 

1 1 

V~2T 45\/l4 ~ 105 y/2 315~/ 



Z1/ n 



(All) 



(A.12) 



As for the symmetric matrices, their nonzero elements are 

/ 



-,[000] 



H? 



26 
35 



22 



1 



105 3V10 
22 8 1 



45 



315V22 



V 105 105 7V10 45 Vu 315\/2 315\/22 / 



-.[on] 



f 3 - -lo 

5 10 
1 4 



\ 



5V14 
1 1 



V 10 15 3^ 5VU/ 



-,[022] 



(- 

2 3 
V 2 



-foo 





(A.13a) 



(A.13b) 
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Table A.2 

Symmetry and spectral leakage I (4.26) of the matrices j 1 ^!^] an( j 



rp[kd\d<2] 

H? 





_,[fed 1 d 2 ] 
1 Hi 






~,[fed 1 d 2 ] 




\kd\d2\ 


Symmetry 


I 


[fcdicfe] 


Symmetry 


I 





S 


4 





S 


4 


10 


s 


5 


10 


N/A 


3 


2 


S 


6 


Oil 


S 


2 


Oil 


s 


2 


2 2 


s 











10 


s 


5 








110 


N/A 


4 








111 


S 


3 








2 


S 


6 








2 11 


s 


4 



-,[200] 



-,[100] 



94 
315 



26 
35 



22 



V 105 105 7\/l0 45 \/l4 315^2 315^22 / 



.[in] 

H? 



16 



105 

8 


3n 


105 




/ 


1 







To 


1 


2 


V 10 


15 



45 



315V22 



1 



1 



14 



315 21V10 
16 4 1 



495 



V2 
315 



1 



35 \/2 
1 

15 35v^/ 

7= 

4095V22 

1 2 



2 J -2- \ 

15 1 



9009 



2 o/2 
13 Z V 15 



V 315 315 63VT0 165 s/li 693\/2 1365^22 3465 9009 / 

.ii i 





,[211] 
H? 



3 1 
35 70 







2 \ 

V ii 4 



1 16 



15V14 
1 



V2 



105 
V ii 



V 70 105 7VT0 15 VU 105 105 / 
again for to < 2 and, in each case, n < I. 



(A.13c) 



(A.13d) 



(A.13e) 



(A.13f) 



A.3. The Matrices T^ 1 and T^ffl 



Regarding the matrices T^f 1 ^ € 



and (4.12) for m > 4, leads to 



N b xN u 
-,[fcdid 2 ] 



(4.31), the relation pb m — DA^J_3, which follows from (4.10c) 



-,[fe(di + l)d 2 ] 
Hi 



. Therefore, given the results in §A.l, 



- m— 3,n 

we only need to evaluate the elements in rows 1-3. In the main text we make use of the matrices with 
[kdKh] = [000], [001], [011], [012], [100], and [111]. Among these matrices, T l °\% has no nonzero elements 



in its first three rows, and the only corresponding nonzero element of T^i^ 2 is 
Moreover, we have 



-,[011] 

mm 



-lb- 1 ' 2 . 



13 
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-,[000] 

mm 



3\/l0 15\/l4 



3\/l0 15\/Ii 







7,(100] 

1 H^m 



l l 



105\/3 / 
V2 „ 



-,[001] 

h 1 ^ 



3V10 
1 



3V10 



(A.14a) 



5\/21 / 



21-s/lO 15\/l4 315 
1 1 V2 



21 a/10 15\/14 315 
^= 



^[iii] 



3\/10 



3V10 



(A.14b) 



5\/21/ 



15V21 105\/33/ 

As for the T^ffl matrices (4.31), one can deduce from (4.15) the relation \T l ^ 
where n > 3. Thus, it suffices to write down the nonzero elements in their first two columns, namely 

/ 3 7 1 3 





rp[kd\d-2] 


- ran 





^[000] 

± mm 



7,(001] 

± mm 



(I 1_V6 Q 

2 2 5 
111 1 



10 10 V6 7VT0 

2 1 Vs V 5 
V 15 5 5 21 

\ T 



1 

105\/2 
1 



15\/14 105 \/2/ 



1 



5\/!I 
1 



7,(011] 
1 mm 



V 6 6 5\/6 3\/l0 5\/l4/ 



o o -L ' 



10 



7,(100] 



/in 

30 30 



1 



\/6 a/10 

- \ T 
2 J 4- \ 



35 





15 35 



3VT0 45 
1 



7,(111] 



2 

5 

21~ 

T 



315 

2,/^T 
V ii 



0^—0 

5 



1 1 



5V14 
3 



45\/14 105 \/2 315 / 
/ 



7,(012] 
- 1 if iff? 







V 6 6 5 3 5V14 / 
A. 4. The S and C Matrices 



1 1 [3 
V~2 2 V2 / 



(A.15a) 



(A.15b) 



(A.15c) 



(A.15d) 



Consider first the N x N real matrices S'Lr 1 2 ' and cL/ 2 ', where 



»[ddid 2 ] 
& fj-r 



= ((D^)D d2 AM,D dl AH) 0i ^ 



~ [dditZ 2 



= ((D d Cffs )D <i2 AW,D <il AW) ,, } . (A.16) 



Since, as can be checked from (4.10), the polynomial degree of A^ is p = N + 2r — 1 and (4.38) holds for 
polynomial integrands of degree 2G — 1, it follows that 



G > \(2p + 1 - di - d 2 )/2] 



(A.17) 



is sufficient to evaluate (A.16) exactly using Mach's quadrature scheme (4.38) (see Remark 16 below). 
Specifically, introducing the differentiation matrices A [d] e R GxN , where \A [d] ] = D A^cffc 1 ), thc 
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diagonal weight matrix p e M GxG with [p]kk = P^ck > an( ^ m& king use of the symmetry property Vm(— £) 
(— l) m+1 Vm (£), leads to the expressions 



a [ddid 2 



~[dd 1 d 2 ] 

(_/ rrr 



^ _ ^_]^m+n+d+di+d 2 

2 F « 



l + (-l) 



(A.18a) 
(A.18b) 



In order to evaluate the corresponding matrices for the and Hf(f2) bases, we require, in addition 

to A^ (in these cases defined in terms of the fi m and v n polynomials), the differentiation matrices A^ <G 
R GxN , given by 



• d 



DMn(-Cfe), H\Q) basis, 
lD^„(-CS ] ), basis, 
as the nodal shape functions do not have definite symmetry about £ = 0. Note that the degree of /zat and 
j^v is now p = A — 1 and p = N + 1, respectively (see Propositions 2 and 3), and the quadrature order 

G (A. 17) must be modified accordingly. Introducing S H i and C H i , where 

[e/did 2 ] 



~[ddid 2 ] 



= ((I) d S H )I) d2 Pn,~D dl p m ) j 2 , 



' H 1 



= ((D c ff£ )D V«,D Vm) 0/ }, (A.20) 



we obtain 



5^ ld21 = ^A^fpA^ - (-l) d (^ [dll ) T p^ [<i21 ) if|/2, 
d [dd ld2] = ^ A [ dl] y pA [ d2] + { _ 1)d{A ldl] fpA ld2] ^ H*/2, 



(A.21a) 
(A.21b) 



[dd 1 d 2 



and analogous expressions for S ff 2 and C H 2 . We remark that relations similar to (4.24) also apply 
for the matrices in (A.20), and can be used to economize on computational and coding effort. Taking into 
account (4.37), the matrices defined in (4.39a) and (4.39b) follow from 

( G [dd 1 d 2 



X 



smh(H z z ) 



y[ddid 2 



H 2 

[ddid 2 
H? 



channel problems, 
film problems, 



X a[ddid 2 
sinh(H z z ) bb 



(A.22) 



where the matrix dimensions are respectively set to N u x N u and Nf, x Nf,, and the quadrature order 
Gsatisfies (A.17) for the given N u and N b (see Tabic 1). The matrices C [ddld2] and C [dd ld2] can be obtained 
in a similar manner. 

The Nf, x N u matrices g^d ld2 ] in (4.39c), and the corresponding C b ddld2 \ are evaluated by means of a 
small modification of the method described above. Specifically, setting G > \(p u +pb + 1 — d\ — ^2)/2] , where 
p u and pb are respectively the polynomial degrees of the velocity and magnetic-field bases, we compute the 
G x N u differentiation matrices 



A [d] 



D 2 AH channel prob. 

D d2 ^(Cfe)> nlm prob. 



D d2 APl(-£if| J ), channel prob. 



kn 



^Vn{-ti [ G,h film prob. 



and the G x Nh matrices 



.[d] 



kn 



kn 



Then, using (4.37), we obtain 

[ddid 2 ] _ cosh(H z z ) 



bu 



^bu 



{{A^Y pA^- { -lf(Aty pA^ 



H z (cosh(H z ) - 1) 

^Hh z z ) ( s [dl] y [d2] d / 

^(cosh(^)-l) / P " +( 1J I 



^2] 



(A.23) 
(A.24) 

(A.25) 
(A.26) 
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Remark 16. Mach's algorithm [34] for the a n and b n coefficients (with n € {0, 1, 2, ... , G + l}) of polynomi- 
als orthogonal with respect to the weight function e 11 ^ consists of two parts. For n < Go := mm{[H^\, G+l} 
the coefficients are evaluated algebraically, while if G + 1 > [H^] an iterative procedure is used for n > Go- 
We observed that for the typical and G used in our linear-stability schemes (both of which are signifi- 
cantly larger than the ones considered in Mach's paper), the quadrature knots £ G £ and weights p G | , which 
follow from the eigenvalues and eigenvectors of the Jacobian matrix J constructed from {o n } and {b n } [51], 
are more accurately computed if the iterative procedure is employed for all n. Moreover, using a specialized 
solver for symmetric tridiagonal matrices (e.g. the LAPACK routine DSTEV [48]), rather than a generic 
one, enhances the stability of the computation at large G. Regarding the algorithm's large- behavior, in 
64-bit arithmetic the weight calculation overflows at around — 700. This limitation can be mitigated by 
increasing the arithmetic precision, but doing so is significantly more complicated than in the case of the 
LGL method (see Remark 12), as it involves porting the routines for the J eigenproblem. 

Appendix B. Eigenvalues of Selected Film and Channel Problems 

This appendix contains tables of eigenvalues for the stability problems studied in §5.1. In each case, 
the eigenproblem (4.3) has been solved using the QZ algorithm, and the resulting complex phase velocity 
c — i"f/a is listed in order of decreasing Im(c). In the examples where the spectrum exhibits the A, P, and 
S branches (Tables B.1-B.5 and B.8) the modes are also labeled in order of decreasing Im(c) within their 
respective families. In Tables B.6 and B.7, hydrodynamic and magnetic modes are respectively labeled H 
and M. In addition to c, Tables B.2 and B.6-B.8 also display the modal energies (2.15). All problems with 
H z > (Tables B.3, B.4 and B.6-B.8) have the Hartmann profiles (2.12). In these cases, the stiffness matrix 
K has been computed by means of the exact-quadrature method (Eqs. (4.40) and (4.41)), aside from the 
inductionless problems in Table B.3, where LGL quadrature (4.42a) has also been used. The free-surface 
parameters are Pg = 1.10 x 10~ 4 and Oh = 3.14 x 10~ 4 for all film problems. 
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Table B.l 

Complex phase velocity of the 33 least stable modes of non- MHD channel flow at Re = 10 4 , a = 1, and p u = 500. E and O 



respectively denote 


even 


and odd modes. The underlined digits 


differ from Table VII in Kirchncr [6]. 






Symmetry 


c 


1 


Ai 


E 


2.375264888204708E - 01 + 3.739670622977800E - 03 i 


2 


Pi 


O 


9.646309154505980E - 01 - 3.516727763102788E - 02 i 


3 


P 2 


E 


9.646425100392813E - 01 - 3.518658379244306E - 02 i 


4 


A 2 


O 


2.772043438088044E - 01 - 5.089872725696847E - 02 i 


5 


Ps 


O 


9.363165358813226E - 01 - 6.320149583992000E - 02 i 


6 


P 4 


E 


9.363517811647262E - 01 - 6.325156907426749E - 02 i 


7 


Ps 


O 


9.079830546294242E - 01 - 9.122273543365197E - 02 i 


8 


Pe 


E 


9. 080563344920716 E - 01 - 9.131286177904398E - 02 i 


9 


P 7 


O 


8.796272922073755E - 01 - 1.192328526197428E - 01 i 


10 


P 8 


E 


8.797556958146369E - 01 - 1.193707310085290E - 01 i 


11 


A 3 


E 


3.491068201236165E - 01 - 1.245019775533875E - 01 i 


12 


A 4 


O 


4.163510155757348E - 01 - 1. 382265253008630 E - 01 i 


13 


Pg 


O 


8.512458401242534E - 01 - 1.472339290757578E - 01 i 


14 


Pio 


E 


8.514493818793474E - 01 - 1.474256007531364E - 01 i 


15 


Pn 


O 


8.228350406948775E - 01 - 1.752286786602681E - 01 i 


16 


Pl2 


E 


8.231369612662293E - 01 - 1.754780735526174E - 01 i 


17 


A 5 


E 


1.900592493682310E - 01 - 1.828219254122344E - 01 i 


18 


A 6 


O 


2.127257823532073E - 01 - 1.993606947537197E - 01 i 


19 


Pl3 


O 


7.943883849443799E - 01 - 2.032206650247992E - 01 i 


20 


Pl4 


E 


7.948183878257583E - 01 - 2.035291440392746E - 01 i 


21 


A 7 


O 


5.320452087682050E - 01 - 2.064652191000982E - 01 i 


22 


A s 


E 


4.749011869521779E - 01 - 2.087312200487454E - 01 i 


23 


Pis 


O 


7.658768104770047E - 01 - 2.311859867813260E - 01 i 


24 


Pl6 


E 


7.664940762955391E - 01 - 2.315850738470669E - 01 i 


25 


Ag 


E 


3.684984783493122E - 01 - 2.388248317187859E - 01 i 


26 


Pl7 


O 


7.374157634158677E - 01 - 2.587170762850995E - 01 i 


27 


Pl8 


E 


7.381150135550412E - 01 - 2.596918833894374E - 01 i 


28 


Aio 


O 


6.367193719487207E - 01 - 2.598857151068238E - 01 i 


29 


An 


O 


3.839876109047478E - 01 - 2.651064996075768E - 01 i 


30 


A12 


E 


5.872129329806185E - 01 - 2.671617095882041E - 01 i 


31 


Pl9 


O 


7.123158603613657E - 01 - 2.855147362764390E - 01 i 


32 


A13 


E 


5.129162044858087E - 01 - 2.866250415809010E - 01 i 


33 


P20 


E 


7.088746527386480E - 01 - 2.876553928005440E - 01 i 
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Table B.2 

Complex phase velocity and free-surface energy of the 25 least stable modes of non-MHD film problems at a = 1, p u = 500, 

and fie = 10 4 , 3 X 10 4 

c E a /E 

Re = 10 4 



1 


F 


1.679075867077238E 


+ 


00 


— 


5.309992786818418E 


— 


04 i 


5.91709761E 


— 


01 


2 


U/Ai 


-1.643176086760401E 


- 


01 


- 


1.185929268008322E 


- 


02 i 


3.23193182E 


- 


01 


3 


Pi 


9.645730769959389E 


— 


01 


— 


3.500080786980365E 


— 


02 i 


7.72375828E 


— 


03 


4 


P 2 


9.361509487313044E 


— 


01 


— 


6.277440588217570E 


— 


02 i 


7.42572688E 


— 


03 


5 


A 2 


2.741497879597380E 


— 


01 


— 


6.654001933197072E 


— 


02 i 


1.87107321E 


— 


01 


6 


Pa 


9.07768479757351 IE 


- 


01 


- 


9.050782178609683E 


— 


02 i 


4.22349138E 


- 


03 


7 


P 4 


8.792256199009758E 


— 


01 


— 


1.181138603404709E 


— 


01 i 


1.74071971E 


— 


03 


8 


Ps 


8.508843382551756E 




01 


— 


1.457778119522005E 




01 i 


5.69264153E 




04 


9 




4.227712274179599E 




01 




1. 512473936731 166E 




01 i 


4.81242748E 




02 


10 


P 6 


8.221887607323224E 


- 


01 


- 


1.732165968799022E 


- 


01 i 


1.66249755E 


- 


04 


11 


A 4 


2.225816434327771E 


— 


01 


— 


1.922767991882933E 


— 


01 i 


2.02884537E 


— 


02 


12 


P 7 


7.939378206791210E 


— 


01 




2.008584356742881E 


— 


01 i 


4.29839279E 


— 


05 


13 


A 5 


5.422180389789260E 


— 


01 




2.120033027735726E 


— 


01 i 


1.39410132E 


— 


03 


14 


Ps 


7.650320024922520E 


- 


01 




2.280604498178582E 


- 


01 i 


1.09648160E 


- 


05 


15 


P 9 


7.371670485488192E 


— 


01 




2.552883380859681E 


— 


01 i 


2.58781548E 


— 


06 


16 


A 6 


6.456425645536397E 


— 


01 




2.592036510454908E 


— 


01 i 


1.42677517E 


— 


05 


17 


A 7 


3.909888051449190E 


— 


01 




2.598302952800038E 


— 


01 i 


3.52648568E 


— 


04 


18 


Pio 


7.119981885525993E 


— 


01 




2.82416251 1392642E 


- 


01 i 


5.54432313E 


- 


07 


19 


A 8 


5.269953093034085E 


— 


01 




3.098257193477898E 


— 


01 i 


2.79614942E 


— 


06 


20 


Pn 


6.932307437171046E 


- 


01 




3.181523206709904E 


- 


01 i 


4.79906211E 


- 


08 


21 


A 9 


6.454839213372529E 


— 


01 




3.482488429410512E 


— 


01 i 


1.35560506E 


— 


08 


22 


Si 


6.772395513979280E 


— 


01 




3.643694475218588E 


— 


01 i 


1.46995919E 


— 


09 


23 


s 2 


6.739703032625379E 


— 


01 




4.124910427107358E 


— 


01 i 


2.36304963E 


— 


11 


24 


s 3 


6.728003310661427E 


- 


01 




4.593456697575948E 


- 


01 i 


3.24310075E 


— 


12 


25 


S 4 


6.719889441665264E 


- 


01 




5.076462584341503E 


- 


01 i 


5.83055870E 


— 


12 


Re = 3 X 10 4 
























1 


Ai 


1.753678323257369E 


— 


01 


+ 


7.984943826436932E 


— 


03 i 


7.47974380E 


— 


02 


2 


F 


1.177455429168835E 


+ 


00 


+ 5.244714510238900E 


— 


05 i 


6.57040491E 


— 


01 


3 


Pi 


9.794482099282948E 


— 


01 




1.996753274816041E 


— 


02 i 


2.93445909E 


— 


02 


4 


P 2 


9.628898245015168E 


— 


01 




3.562268032337058E 


— 


02 i 


2.68859469E 


— 


02 


5 


Pa 


9.463998634683490E 


— 


01 




5.116697097189563E 


— 


02 i 


1.48771717E 


— 


02 


6 


P 4 


9.297557437185264E 


— 


01 




6.645791488698941E 


— 


02 i 


6.13372852E 


— 


03 


7 


Ps 


9.134687978618840E 


— 


01 




8.177074836942153E 


— 


02 i 


2.00134460E 


— 


03 


8 


A 2 


2.440029358172409E 


— 


01 




8.299437421776400E 


— 


02 i 


7.33186450E 


— 


02 


9 


Pe 


8.969522064288651E 




01 




9.663007449652813E 




02 i 


6.10937302E 




04 


10 


P 7 


8.812525388014069E 




01 




1.119021889404179E 




01 i 


1.56175054E 




04 


11 


Ps 


8.652977734764966E 




01 




1.266004433067371E 




01 i 


4.11484076E 




05 


12 


A 3 


1.319105765324801E 




01 




1.329667089050516E 




01 i 


1.68098585E 




02 


13 


P 9 


8.501964700857818E 




01 




1.425578808374347E 




01 i 


8.41602738E 




06 


14 


A 4 


3.321812125705384E 




01 




1.468796725488414E 




01 i 


7.62560468E 




03 


15 


Pio 


8.345314564493534E 




01 




1.577157358618176E 




01 i 


1.90835305E 




06 


16 


Pn 


8.192521379852682E 




01 




1.744515939429078E 




01 i 


3.31647400E 




07 


17 


A 5 


2.567572382646105E 




01 




1.759022308300596E 




01 i 


1.26817844E 




03 


18 


Pl2 


8.033985922501681E 




01 




1.898886430093124E 




01 i 


6.96533353E 




08 


19 


A 6 


4.122424067140122E 




01 




1.918160189923045E 




01 i 


1.19508181E 




04 


20 


Pia 


7.877723662547554E 




01 




2.069276983777561E 




01 i 


1.13100490E 




08 


21 


A 7 


3.600810105077645E 




01 




2.157286814750727E 




01 i 


2.06150843E 




05 


22 


Pl4 


7.717299604876217E 




01 




2.224013373973636E 




01 i 


2.31395673E 




09 


23 


A 8 


4.858154352656867E 




01 




2.273008596491542E 




01 i 


1.45547586E 




06 


24 


Pl5 


7.558940381686612E 




01 




2.395471826819909E 




01 i 


3.63694162E 




10 


25 


A 9 


4.483456535355433E 




01 




2.515881911262335E 




01 i 


1.80327316E 




07 
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Table B.3 



Complex phase velocity of the 25 least stable modes of inductionless film problems at Re = 3 X 10 4 , H x = 0, a = 1, and 
p u = 500 and H z = 14, 100. The mean basic velocities (2.13) arc (U) = 0.92857 (H z = 14) and (U) = 0.99000 (H z = 100) 

c (exact) c (LGL) 





= 14 




















1 


F 


1.264889076407188E 


+ 00 


- 2.709480082539037E 


- 03 i 


1.264889076407183E 


+ 00 


- 


2.709480082550979E 


- 03 i 


2 


Pi 


9.992558947033491E 


- 01 


- 8.278625506899678E 


- 03 i 


9.992558947033513E 


- 01 


- 


8.278625506898540E 


- 03 i 


3 


P 2 


9.983412661330004E 


-01 


- 1.153953543381368E 


- 02 i 


9.983412661330017E 


- 01 


- 


1.153953543381040E 


- 02 i 


4 


P 3 


9.971872338663629E 


- 01 


- 1.574118837369919E 


- 02 i 


9.971872338663625E 


- 01 


— 


1.574118837369889E 


- 02 i 


5 


P 4 


9.958459023889774E 


- 01 


- 2.082236225171166E 


- 02 i 


9.958459023889762E 


- 01 


— 


2.082236225171246E 


- 02 i 


6 


Ai 


7.299947887448923E 


- 01 


- 2.612159002726035E 


- 02 i 


7.299947887448834E 


- 01 


— 


2.612159002730650E 


- 02 i 


7 


P 5 


9.943367200567261E 


- 01 


- 2.672287190234764E 


- 02 i 


9.943367200567290E 


- 01 


- 


2.672287190234856E 


- 02 i 


8 


P 6 


9.926693788193111E 


- 01 


- 3.340152602696682E 


- 02 i 


9.926693788193063E 


- 01 


- 


3.340152602697706E 


- 02 i 


9 


P 7 


9.908542144571485E 


- 01 


- 4.082332709314562E 


- 02 i 


9.908542144571193E 


- 01 


— 


4.082332709314013E 


- 02 i 


10 


P 8 


9.888944056193677E 


- 01 


- 4.896621663250805E 


- 02 i 


9.888944056193839E 


- 01 




4.896621663247380E 


- 02 i 


11 


Pg 


9.867999867265210E 


- 01 


- 5.780362683832217E 


- 02 i 


9.867999867264968E 


- 01 




5.780362683834655E 


- 02 i 


12 


Pio 


9.845690593545204E 


-01 


- 6.732496184323691E 


- 02 i 


9.845690593544832E 


-01 




6.732496184304332E 


- 02 i 


13 


Pu 


9.822142176490877E 


-01 


- 7.750474943687088E 


- 02 i 


9.822142176494953E 


- 01 




7.750474943687220E 


- 02 i 


14 


Pl2 


9.797277611242884E 


- 01 


- 8.834259012386456E 


- 02 i 


9.797277611242368E 


- 01 




8.834259012432207E 


- 02 i 


15 


Pl3 


9.771270171142551E 


- 01 


- 9.980938378195403E 


- 02 i 


9.771270171138522E 


- 01 




9.980938378193582E 


- 02 i 


16 


A 2 


2.359311079687786E 


- 01 


- 1.105756921787289E 


- 01 i 


2.359311079687282E 


- 01 




1.105756921787712E 


- 01 i 


17 


Pl4 


9.743971849585732E 


- 01 


- 1.119159843147402E 


- 01 i 


9.743971849586051E 


- 01 




1.119159843143422E 


- 01 i 


18 


PlB 


9.715618781861247E 


- 01 


- 1. 24625571 1002706E 


- 01 i 


9.715618781862945E 


- 01 




1. 24625571 1003195E 


- 01 i 


19 


Pl6 


9.685976062545173E 


- 01 


- 1.379622724905363E 


- 01 i 


9.685976062546490E 


- 01 




1.379622724893477E 


- 01 i 


20 


Pit 


9.655345585233294E 


- 01 


- 1.518775344849239E 


- 01 i 


9.655345585280525E 


- 01 




1.518775344852824E 


- 01 i 


21 


Pis 


9.623381642628556E 


-01 


- 1.664108150280590E 


- 01 i 


9.623381642634273E 


-01 




1.664108150381094E 


-Oli 


22 


Pl9 


9.590419783471862E 


- 01 


- 1.814972823595940E 


- 01 i 


9.590419783333189E 


- 01 




1.814972823657535E 


- Oli 


23 


P20 


9.555914363432764E 


- 01 


- 1.971912711398718E 


- 01 i 


9.555914363242426E 


- 01 




1.971912711279786E 


- Oli 


24 


P21 


9.520009797772009E 


- 01 


- 2.133960841540300E 


- 01 i 


9.520009797733312E 


- 01 




2.133960841184933E 


- Oli 


25 


P22 


9.481381471570949E 


- 01 


- 2.301291550139524E 


- 01 i 


9.481381472091938E 


- 01 




2.301291549732140E 


- Oli 



H z 


= 100 


























1 


F 


1.244657312459241E 


+ 


00 




1.276484109616633E 




Oli 


1.244657312459336E 


+ 00 




1.276484109616758E 


- Oli 


2 


Ai 


7.527712613119597E 




01 




1.309887753320078E 




Oli 


7.527712613118247E 


- 01 




1.309887753319950E 


- Oli 


3 


Pi 


9.992315130288280E 




01 




3.101349562479781E 




Oli 


9.992315130288494E 


- 01 




3.101349562479994E 


-Oli 


4 


P 2 


9.997105114818259E 




01 




3.284913791218512E 




Oli 


9.997105114818199E 


- 01 




3.284913791218546E 


- Oli 


5 


P 3 


9.997646389857329E 




01 




3.342041754324304E 




Oli 


9.997646389857272E 


- 01 




3.342041754324196E 


- Oli 


6 


P 4 


9.997173589263894E 




01 




3.383276198992341E 




Oli 


9.997173589263828E 


-01 




3.383276198992299E 


- Oli 


7 


Ps 


9.996252105204291E 




01 




3.424685385378852E 




Oli 


9.996252105204261E 


- 01 




3.424685385378918E 


- Oli 


8 


Pe 


9.995023725166139E 




01 




3.470293583958746E 




Oli 


9.995023725166103E 


- 01 




3.470293583958675E 


- Oli 


9 


P 7 


9.993543622637743E 




01 




3.521441474868715E 




Oli 


9.993543622637748E 


- 01 




3.521441474868756E 


- Oli 


10 


P 8 


9.991827574316863E 




01 




3.578762247930776E 




Oli 


9.991827574316822E 


- 01 




3.578762247930801E 


- Oli 


11 


P9 


9.989895462677130E 




01 




3.642452387564944E 




Oli 


9.989895462677091E 


- 01 




3.642452387564947E 


- Oli 


12 


Pio 


9.987747691956044E 




01 




3.712741930516902E 




Oli 


9.987747691956066E 


- 01 




3.712741930516907E 


- Oli 


13 


P11 


9.985398440924721E 




01 




3.789594000704576E 




Oli 


9.985398440924742E 


- 01 




3.789594000704578E 


- Oli 


14 


P12 


9.982845714861840E 




01 




3.873167576511015E 




Oli 


9.982845714861908E 


- 01 




3.873167576510981E 


- Oli 


15 


Pl3 


9.980101265689799E 




01 




3.963348225337110E 




Oli 


9.980101265689778E 


- 01 




3.963348225337082E 


- Oli 


16 


Pl4 


9.977164108934937E 




01 




4.060287849173049E 




Oli 


9.977164108935060E 


- 01 




4.060287849173056E 


- Oli 


17 


PlB 


9.974043604226334E 




01 




4.163838704651704E 




Oli 


9.974043604226298E 


- 01 




4.163838704651684E 


- Oli 


18 


Pl6 


9.970741485284107E 




01 




4.274161109594051E 




Oli 


9.970741485284181E 


- 01 




4.274161109594042E 


- Oli 


19 


Pit 


9.967264025347967E 




01 




4.391091672124338E 




Oli 


9.967264025347992E 


- 01 




4.391091672124348E 


- Oli 


20 


Pis 


9.963617205307972E 




01 




4.514802977904322E 




Oli 


9.963617205308041E 


- 01 




4.514802977904337E 


- Oli 


21 


Pl9 


9.959803285139238E 




01 




4.645125250143292E 




Oli 


9.959803285139230E 


- 01 




4.645125250143295E 


- Oli 


22 


P20 


9.955834489583942E 




01 




4.782243941461477E 




Oli 


9.955834489584048E 


- 01 




4.782243941461505E 


- Oli 


23 


P21 


9.951708277082375E 




01 




4.925988988424788E 




Oli 


9.951708277082342E 


- 01 




4.925988988424814E 


- Oli 


24 


P22 


9.947446209812095E 




01 




5.076558185209421E 




Oli 


9.947446209812166E 


- 01 




5.076558185209454E 


- Oli 


25 


P23 


9.943040892241203E 




01 




5.233785783714180E 




Oli 


9.943040892241175E 


- 01 




5.233785783714235E 


- Oli 
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Table B.4 



Complex phase velocity 
tf z /tan(l°) = 802.06, a 


c of the 33 least stable modes of an 
= 1, and p u = 500 


inductionless channel problem at 


Re = 10 4 , H z = 14, H x = 


1 


Pi 


9.993700195546175E 


- 01 


- 9.721004530793138E - 02 i 


2 


P 2 


9.976612096433158E 


- 01 


- 1.030917967244671E - 01 i 


3 


P 3 


9.951774765278338E 


- 01 


- 1.123185821184923E - 01 i 


4 


P 4 


9.921089021026304E 


- 01 


- 1.244427518042686E - 01 i 


5 


P 5 


9.885576846065738E 


- 01 


- 1.391394166517699E - 01 i 


6 


P 6 


9.845858924007105E 


- 01 


- 1.561699942143518E - 01 i 


7 


Pt 


9.802375881163239E 


- 01 


- 1.753455921597468E - 01 i 


8 


Ps 


9.755530144308208E 


- 01 


- 1.965030475951849E - 01 i 


9 


Pg 


9.705872940874839E 


- 01 


- 2.194897048867427E - 01 i 


10 


Pio 


9.654406146046013E 


- 01 


- 2.441635533378584E - 01 i 


11 


Pn 


9.994941998893740E 


- 01 


- 2.565041224611593E - 01 i 


12 


Pl2 


9.981140152080975E 


- 01 


- 2.611657657075488E - 01 i 


13 


Ai 


6.455703277951147E 


-01 


- 2.663992874104253E - 01 i 


14 


Pl3 


9.960868515182631E 


-01 


- 2.685202959225929E - 01 i 


15 


Pl4 


9.602887830494103E 


-01 


- 2.704320708227363E - 01 i 


16 


Pis 


9.935265121923997E 


- 01 


- 2.782735041820532E - 01 i 


17 


Pl6 


9.905224108407159E 


- 01 


- 2.902400425540363E - 01 i 


18 


Pl7 


9.553532017423860E 


- 01 


- 2.983346408422508E - 01 i 


19 


Pis 


9.870642174734445E 


- 01 


- 3.042371579569956E - 01 i 


20 


Pl9 


9.831554449489976E 


- 01 


- 3.203013628990474E - 01 i 


21 


P20 


9.507606261162745E 


- 01 


- 3.280616396399530E - 01 i 


22 


P21 


9.789956131371971E 


-01 


- 3.382151445827348E - 01 i 


23 


P22 


9.742217743785576E 


-01 


- 3.577645640279395E - 01 i 


24 


P 23 


9.464376866330062E 


- 01 


- 3.598000886532416E - 01 i 


25 


P 2 4 


9.689648687655241E 


- 01 


- 3.794158635141329E - 01 i 


26 


P25 


9.421826293777356E 


- 01 


- 3.935372132859687E - 01 i 


27 


P26 


9.640327742098715E 


- 01 


- 4.025289154041328E - 01 i 


28 


P27 


9.587353922299208E 


- 01 


- 4.262057624377996E - 01 i 


29 


P28 


9.379722615037215E 


- 01 


- 4.289398118832719E - 01 i 


30 


P 29 


9.528921401648801E 


- 01 


- 4.518940856671111E - 01 i 


31 


A 2 


5.997926751489199E 


- 01 


- 4.602178850424406E - 01 i 


32 


P30 


9.343088725666473E 


-01 


- 4.655346048131246E - 01 i 


33 


P31 


9.488034208886676E 


- 01 


- 4.805405779007103E - 01 i 



[30] L. N. Trefethen, M. Embree, Spectra and Pscudospcctra: The Behavior of Nonnormal Matrices and Operators, Princeton 

University Press, Princeton, 2005. 
[31] P. J. Schmid, ct al., A study of eigenvalue sensitivity for hydrodynamic stability operators, Theoret. Comput. Fluid 

Dynamics 4 (1993) 227. 

[32] L. M. Mack, A numerical study of the temporal eigenvalue spectrum of the Blasius boundary layer, J. Fluid Mcch. 73 (3) 
(1976) 497. 

[33] S. C. Rcddy, P. J. Schmid, D. S. Hcnningson, Pscudospcctra of the Orr-Sommcrfcld operator, SIAM J. Appl. Math 53 (1) 
(1993) 15. 

[34] R. Mach, Orthogonal polynomials with exponential weight in a finite interval and application to the optical model, J. 

Math. Phys. 25 (7) (1984) 2186. 
[35] P. G. Ciarlet, Basic error estimates for elliptic problems, in: P. G. Ciarlct, J. L. Lions (Eds.), Finite Element Methods 

(Part 1 ), 1st Edition, Vol. II of Handbook of Numerical Analysis, Elsevier Science B.V., Amsterdam, 1991, p. 17. 
[36] M. O. Deville, P. F. Fischer, E. H. Mund, High-Order Methods for Incompressible Fluid Flow, Vol. 9 of Cambridge 

Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2002. 
[37] U. Banerjeee, J. E. Osborn, Estimation of the effect of numerical integration in finite element eigenvalue approximation, 

Numer. Math. 56 (1990) 735. 

[38] I. Babuska, J. Osborn, Eigenvalue problems, in: P. G. Ciarlet, J. L. Lions (Eds.), Finite Element Methods (Part 1 ), 1st 

Edition, Vol. II of Handbook of Numerical Analysis, Elsevier Science B.V., Amsterdam, 1991, p. 641. 
[39] J. A. Shcrcliff, A Textbook of Magnctohydrodynamics, Pcrgamon Press, Oxford, 1965. 

[40] J. T. Stuart, On the stability of viscous flow between parallel planes in the presence of a co-planar magnetic field, Proc. 

Roy. Soc. A. 221 (1145) (1954) 189. 
[41] G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967. 
[42] A. Hof, et al., Onset of oscillatory convection in molten gallium, J. Fluid Mech. 515 (2004) 391. 
[43] S. C. Hardy, The surface tension of liquid gallium, J. Cryst. Growth 71 (1985) 602. 

[44] V. Kolevzon, Anomalous temperature dependence of the surface tension and capillary waves at a liquid gallium surface, 
J. Phys.: Condcns. Matter 11 (1999) 8785. 



53 



Table B.5 



Complex phase velocity c of the 25 least stable magnetic modes of film MHD flow with zero background magnetic field 
(H x = H z = 0) at Re = 10 4 , Pm = 1.2, a = 1, and p u = Pb = 500. The numbering in the left-hand column takes into account 
the hydrodynamic part of the spectrum, where the eigenvalues are identical to those listed in the Re = 10 4 portion of Table B.2. 



2 


Pmi 


9.931687361486334E 


- 01 - 7.411597582310019E - 03 


4 


Pm 2 


9.675367306486845E 


- 01 - 3.281002191260387E - 02 


6 


Ami 


6.190480900456202E 


- 02 - 3.658789628263018E - 02 


7 


Pm 3 


9.417639706299372E 


- 01 - 5.851812430456833E - 02 


10 


Pm 4 


9.159676389192429E 


- 01 - 8.428165278355801E - 02 


12 


Am 2 


1.920049933116885E 


- 01 - 1.071612998031651E - 01 


13 


Pm 5 


8.901624704645013E 


- 01 - 1.100662088133551E - 01 


15 


Pm 6 


8.643528846481575E 


- 01 - 1.358613270394366E - 01 


18 


Am3 


2.822530721072338E 


- 01 - 1.539360769033924E - 01 


19 


Pm 7 


8.385407224351258E 


- 01 - 1.616626213801015E - 01 


21 


Pm 8 


8.127269041407206E 


- 01 - 1.874678918924152E - 01 


22 


A1114 


3.582077439231074E 


- 01 - 1.914285425627469E - 01 


26 


Pm 9 


7.869119473724944E 


- 01 - 2.132758984008717E - 01 


27 


Am 5 


4.256621499174164E 


- 01 - 2.230812126023304E - 01 
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Table B.6 

Complex phase velocity, free-surface energy and magnetic energy of the 25 least stable modes of film MHD problems at 

Re = 10 4 , Pro = 1.2, H x = 0, a = 1, p u = p b = 500 and H z = 14, 100 
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Table B.7 

Complex phase velocity, free-surface energy, and magnetic energy of the 50 least stable modes of film MHD flow at Re = 10 4 , 
Pm = 1.2, H z = 100, H x = 100/ tan(l°) = 5,729.0, a = 1, and p u = p b = 500 
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48 


H24 


5 


519445729025135E 


- 01 




4 


922805765858015E 




01 i 


1 


18405939E 




09 


5 


00000205E 




01 


49 


M25 


-7 


030941731336009E 


- 01 




5 


149608442491478E 




01 i 


1 


33668718E 




03 


5 


00079798E 




01 


50 


M 26 


2 


689142835219470E 


+ 00 




5 


151046749735411E 




01 i 


1 


93250782E 




09 


4 


99999989E 




01 



56 



Table B.8 

Complex phase velocity, free-surface energy, and magnetic energy of the 25 least stable modes of Pm = 10~ 4 and inductionlcss 



film problems 


at Re = 


10 6 , a = 0.01, H x 




0, 


H 


z ■ 


= 10, and p u = pi, = 


_ i 


500 




























c 












E a /E 








E b /E 




Pm 


= 10~ 4 


































1 


Pi 


9 


827762904552182E 




01 + 1 


379303998591685E 


— 


02 i 


7 


.21784034E 


— 


02 


2 


.15340701E 


- 01 


2 


P 2 


9 


966319464666803E 




01 


— 


1 


729497204355026E 


— 


02 i 


1 


.59127542E 


— 


03 


6 


39646788E 


- 03 


3 


F 


1 


.014164030326762E + 00 


— 


2 


023492335823250E 


— 


02 i 


5 


.55459542E 


— 


02 


2 


.02404403E 


- 01 


4 


M 


8 


963943583270458E 


— 


01 


— 


2 


230611416635310E 


— 


02 i 


2 


.62501396E 


— 


03 


2 


.81923599E 


- 01 


5 


P 3 


9 


.924107060542195E 




01 


— 


2 


791684540382556E 


- 


02 i 


1 


37991443E 


- 


03 


5 


.42646365E 


- 03 


6 


P 4 


9 


.875153580671349E 




01 


— 


4. 


091279696616719E 


— 


02 i 


5 


.85319839E 


— 


04 


2 


.43710490E 


- 03 


7 


P 5 


9 


821430916134596E 


- 


01 




5 


598669702561006E 


- 


02 i 


2 


15341228E 


— 


04 


9 


.67853643E 


- 04 


8 


P 6 


9 


763394322920629E 


— 


01 




7 


303717929168575E 


- 


02 i 


7. 


94688259E 


— 


05 


3 


.87838637E 


- 04 


9 


P 7 


9 


701177299578914E 


— 


01 


- 


9 


198227116387199E 


- 


02 i 


3 


02967004E 


— 


05 


1 


.61289078E 


- 04 


10 


P8 


9 


.634914456358989E 


— 


01 


— 


1 


127436343636547E 


- 


01 i 


1 


.19476180E 


- 


05 


7 


11236630E 


- 05 


11 


P 9 


9 


564737788331985E 


— 


01 


— 


1 


352508362860170E 




01 i 


4 


.84917076E 




06 


3 


.43418962E 


- 05 


12 


Ai 


2 


043753294235331E 


— 


01 


— 


1 


382155401377719E 




01 i 


2 


64374024E 




05 


6 


.76888446E 


- 04 


13 


Pio 


9 


.490769002736642E 


— 


01 


— 


1 


594407787552430E 


- 


01 i 


2 


.01577028E 


- 


06 


1 


.86829537E 


- 05 


14 


Pn 


9 


413059285197984E 


— 


01 


— 


1 


852666320174906E 


- 


01 i 


8 


.55277725E 


- 


07 


1 


.15816101E 


- 05 


15 


Pl2 


9 


.331100318462960E 


— 


01 


— 


2 


126992189563545E 


- 


01 i 


3 


.70391337E 


- 


07 


8 


.07320728E 


- 06 


16 


Pl3 


9 


241628962902693E 


— 


01 


— 


2 


415996137458210E 


— 


01 i 


1 


.66752772E 


— 


07 


6 


.18112296E 


- 06 


17 


A 2 


5 


485856863523154E 


— 


01 


— 


2 


.543247513754485E 


- 


01 i 


2 


.03272277E 


- 


05 


2 


73275033E 


- 04 


18 


Pl4 


9 


.136143986279063E 


— 


01 


— 


2 


701419221344245E 


- 


01 i 


8 


.65851959E 


- 


08 


5 


.13716300E 


- 06 


19 


A 3 


7 


584423091150166E 


— 


01 


— 


2 


936517999528958E 


- 


01 i 


2 


.34972269E 


— 


06 


2 


.65866203E 


- 05 


20 


Si 


9 


.082386758765671E 


— 


01 


— 


2 


949546494082194E 


- 


01 i 


4 


.28420360E 


- 


08 


4 


.41251095E 


- 06 


21 


s 2 


9 


091850939621119E 


— 


01 


— 


3 


278070972139672E 


— 


01 i 


1 


25177303E 


— 


08 


3 


67249686E 


- 06 


22 


s 3 


9 


.085887533412148E 


— 


01 


— 


3 


657495322686370E 


— 


01 i 


3 


52956711E 


— 


09 


3 


.15176871E 


- 06 


23 


s 4 


9 


.075729146748517E 


— 


01 


— 


4. 


061278815870065E 


— 


01 i 


1 


.02704210E 


— 


09 


2 


.76414147E 


- 06 


24 


s 5 


9 


065902759286738E 


— 


01 


— 


4. 


485775139379218E 


- 


01 i 


3 


.16436715E 


- 


10 


2 


.45730082E 


- 06 


25 


s 6 


9 


.057202006968366E 


— 


01 


— 


4. 


929937573773935E 


- 


01 i 


1 


.10787600E 


— 


10 


2 


.20552966E 


- 06 


Zero- Pm 


































1 


F 


1 


003081951615145E 


+ 


00 


— 


3 


558190644980472E 


- 


03 i 


6 


.73619474E 


- 


01 








2 


Pi 


9 


.962565311692007E 


— 


01 


— 


8. 


495142270657967E 


— 


03 i 


2 


86281193E 


— 


01 








3 


P 2 


9 


954753400143420E 


— 


01 


— 


1 


799718315702588E 


- 


02 i 


1 


86586723E 


- 


02 








4 


P 3 


9 


914237146859227E 


— 


01 


— 


2 


854506065914169E 


— 


02 i 


4 


.19130496E 


— 


03 








5 


P 4 


9 


.867952030406444E 


— 


01 


— 


4. 


138105767960570E 


— 


02 i 


1 


11289257E 


— 


03 








6 


P 5 


9 


816310168588849F 




01 


— 


5 


636524362581984E 


— 


02 i 


3 


36531038E 


— 


04 








7 


Pe 


9 


.759749612552397E 




01 


— 


7. 


336441522861509E 


- 


02 i 


1 


.11752029E 


- 


04 








8 


P 7 


9 


698607530025503E 




01 




9 


227263488200604E 


— 


02 i 


3 


.97310832E 


— 


05 








9 


Ps 


9 


.633151 174989847E 




01 




1 


130045908951411E 


— 


01 i 


1 


48716119E 


— 


05 








10 


P 9 


9 


563600657929386E 




01 




1 


354908457302540E 


— 


01 i 


5 


79468653E 


— 


06 








11 


Ai 


2 


040820291301501E 




01 




1 


368283221616525E 


— 


01 i 


2 


.66389530E 


— 


05 








12 


Pio 


9 


490138122621106E 




01 




1 


596757535114703E 




01 i 


2 


33183493E 




06 








13 


Pn 


9 


412830103463788E 




01 




1 


855234968930203E 




01 i 


9 


63548983E 




07 








14 


Pl2 


9 


331097482194961E 




01 




2 


130253760015867E 




01 i 


4 


08075498E 




07 








15 


Pl3 


9 


.241307161531740E 




01 




2 


420975515207121E 




01 i 


1 


80152651E 




07 








16 


A 2 


5 


.4729946561 12545E 




01 




2 


566280626553634E 




01 i 


2 


.03243362E 




05 








17 


Pl4 


9 


130639876833722E 




01 




2 


709666671986019E 




01 i 


9 


.35774907E 




08 








18 


A 3 


7 


599241470438746E 




01 




2 


9250421 13724450E 




01 i 


2 


.52448668E 




06 








19 


Si 


9 


.073118082454555E 




01 




2 


946601335295559E 




01 i 


4 


87403966E 




08 








20 


s 2 


9 


090167819631659E 




01 




3 


275242199968214E 




01 i 


1 


36448737E 




08 








21 


s 3 


9 


085298964747124E 




01 




3 


656433228187246E 




01 i 


3 


77374446E 




09 








22 


s 4 


9 


.075427230252736E 




01 




4. 


060870759697409E 




01 i 


1 


.08690923E 




09 








23 


s 5 


9 


065734350294627E 




01 




4. 


485622185122770E 




01 i 


3 


32509786E 




10 








24 


s 6 


9 


.057105331801750E 




01 




4. 


929883936040092E 




01 i 


1 


.15779448E 




10 








25 


s 7 


9 


049563953966798E 




01 




5 


393404670912730E 




01 i 


5 


11033938E 




11 
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